r  : 

T  7  “ 

''  AD-A119  074  AIR  FORCE  WRIGHT  AERONAUTICAL  LABS  WRIGHT-PATTERSO 

NAVIER-STOKES  SOLUTIONS  FOR  AN  AXlSYMMETRlC  NOZZLE 
MAR  82  G  A  HASEN  , 

UNCLASSIFIED  AFWAL-TR-81-3161 

AFB  OH  F/G  20/4 

IN  A  SUPERSO-- ETC<U 

NL 

- 1 

1 

tm 

II 

II 

II 

II 

II 

1] 

! 

3^ 

;  & 

ffi 

If 

■ _ 

NAVIER-STOKES  SOLUTIONS  FOR  AN  AXISYMMETRIC 
NOZZLE  IN  A  SUPERSONIC  EXTERNAL  STREAM 


Gerald  A.  Hasen,  Captain,  USAF 
Aerodynamics  and  Airframe  Branch 
Aeromechanics  Division 


March  1982 


Final  Report  for  Period  June  1978  -  April  1981 


Approved  for  public  release;  distribution  unlimited 


FLIGHT  DYNAMICS  LABORATORY 
AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 
AIR  FORCE  SYSTEMS  COMMAND 
WRIGHT-PATTERSON  AIR  FORCE  BASE,  OHIO  45433 


NOTICE 


When  Government  drawings ,  specifications,  or  other  data  are  used  for  any  purpose 
other  than  in  connection  with  a  definitely  related  Government  procurement  operation, 
the  United  States  Government  thereby  incurs  no  responsibility  nor  any  obligation 
whatsoever;  and  the  fact  that  the  government  may  have  formulated,  furnisned ,  or  :n 
any  way  supplied  the  said  drawings ,  specifications ,  or  other  data,  is  not  to  oe  re¬ 
garded  by  implication  or  otherwise  as  in  any  manner  licensing  the  holder  or  any 
other  person  or  corporation,  or  conveying  any  rights  or  permission  to  mar.u:  acture 
use,  or  sell  any  patented  invention  that  may  in  any  way  be  related  tr-er-to. 

This  reoort  has  been  reviewed  by  the  Office  of  Public  Arc  airs  f  - ~  0/  mA )  *r.d  .s 
releasable  to  the  National  Technical  Information  Service  (NTIS) .  -->t  NTiS, 

be  available  to  the  general  public,  including  foreign  nations. 


This  technical  report  has  been  reviewed  and  is  approved  r or  pur  * -cation . 


r 


GERALD  A.  HASEN,  Captain,  USAF 
Project  Engineer 


LOWELL  C.  KEEL,  Lt  Col,  USAK 

Chief,  Aerodynamics  and  Airframe  Branch 


FOR  THE  COMMAND. 


Xr.  CHEVALIER,  Colonel,  USAF 
lief.  Aeromechanics  Division 


If  uour  address  has  changed,  if  you  wish  to  be  removed  from  our  mailing  list,  or 
if  the  addressee  is  no  longer  employed  by  you .  organization  please  notify  AiWAL/MMM, 
W-PAFB,  OH  45433  to  help  us  maintain  a  current  mailing  list  . 

Copies  of  this  report  should  not  be  returned  unless  return  is  required  by  security 
considerations ,  contractual  obligations,  or  notice  on  a  specific  document. 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  nat*  Filtered) 

I  REPORT  DOCUMENTATION  PAGE 


t  report  NUMBER 

AFWAL-TR-81-3161 

4  TITLE  (and  'Subtitle) 


2  GOVT  ACCESSION  NO. 

ftO-hulo  7</ 


NAVIER-STOKES  SOLUTIONS  FOR  AN  AX1SYMMETRIC 
NOZZLE  IN  A  SUPERSONIC  EXTERNAL  STREAM 


READ  INSTRUCTIONS 

_ BEFORE  COMPLETING  FORM 

3  RECIPIENT’S  CATALOG  NUMBER 


5  TYPE  OE  REPORT  A  PERIOD  COVERED 

Final  Report 

.Tune  1978  -  April  1981 

6  PERFORMING  ORG.  REPORT  NUMBER 


17  authoro; 


8  CONTRACT  OR  GRANT  NUMBERO; 


Gerald  A.  Hasen,  Captain,  USAF 

’i  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Flight  Dynamics  Laboratory  (AFWAL/FIMM) 

Air  Force  Wright  Aeronautical  Laboratories  (AFSC) 
Wright-Patterson  AFB,  Ohio  45433 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Flight  Dynamics  Laboratory  (AFWAL/FIM) 

Air  Force  Wright  Aeronautical  Laboratories  (AFSC) 
Wright-Patterson  AFB.  Ohio  45433 _ 

14  MONITORING  AGENCY  name  a  ADDRESSfi/  different  from  Controlling  Office) 


10.  PROGRAM  element.  PROJECT.  T  ask 
AREA  A  WORK  UNIT  NUMBERS 

Program  Element  61102F, 
Project  2307,  Task  2307N6, 
Work  Unit  2307N603 _ 

12.  REPORT  DATE 

March  1982 

13  NUMBER  OF  PAGES 

139 

IS.  SECURITY  CLASS,  (ol  this  report) 


16  DISTRIBUTION  STATEMENT  (ol  this  Report) 


Unclassified 

1 5*.  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


Approved  for  public  release;  distribution  unlimited. 


[t7.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Block  20,  if  different  from  Report) 


16  SUPPLEMENTARY  NOTES 


The  material  reported  herein  is  based  on  the  author's  dissertation  submitted 
in  partial  fulfillment  of  the  requirements  for  the  Doctor  of  Philosophy 
degree  at  the  Air  Force  Institute  of  Technology,  Wright-Patterson  AFB,  Ohio. 


19  KEY  WOROS  (Continue  on  reverse  side  if  necessary  and  identify  by  block  number) 

Computational  Fluid  Dynamics  Compressible  Flow 
Numerical  Solution  Turbulent  Flow 

Navier-Stokes  MacCormack's  Explicit  Method 

Axisymmetric  Nozzle  Adaptive  Grid 

Mach  Disc _ 

20  ABSTRACT  (Continue  on  reverse  side  If  necessary  and  Identify  by  block  number) 

^■Numerical  solutions  of  the  Navier-Stokes  equations  are  obtained  for  an  axi¬ 
symmetric  nozzle  in  a' supersonic  external  stream  (M^  =  1.94,  M.i  =  3.0, 

Re<„  =  2.2  xlO6).  Five  jet  pressure  conditions  ranging  from  a  highly  over- 
expanded  case  which  exhibits  a  "Mach  disc  shock  formation  to  a  slightly  under¬ 
expanded  case  are  examined  and  solved  numerically.  MacCormack's  explicit 
method  is  applied  as  the  numerical  algorithm.  An  adaptive  grid  scheme  is 
utilized  in  the  nozzle  wake  to  allow  the  fine  mesh  region  of  the  computational 
grid  to  remain  in  areas  containing  relatively  high  flow  gradients.  Locally\  , 


DD  ,  1473 


EDITION  OF  I  NOV  6S  IS  OBSOLETE 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  rRTicn  Da  la  Emeradl 


_ Uj_  FLED _ 

SECURITY  CL  a  iCATION  OF  This  PAGEflWlan  D.t.  E nfrtdj 


^dependent  eddy  viscosity  modelling  is  applied  in  the  form  of  a  Cebeci-Smith 
two- layer  model  in  the  boundary  layer  regions  on  the  nozzle  walls,  and  a  form 
of  the  Prandtl  mixing  length  model  in  the  nozzle  wake.  A  two-dimensional  wedge 
flat  plate  validation  case  was  computed  using  these  models  with  excellent 
results.  The  computational  solutions  of  the  axisymmetric  nozzle  accurately 
reproduced  the  experimentally  observed  viscous  effects  on  the  nozzle  base 
pressure  and  shock  wave  locations  that  are  caused  by  the  thick  nozzle  base 
annulus.  Correct  transition  was  achieved  numerically  from  regularly  reflected 
shock  waves  at  the  line  of  symmetry  in  the  jet  core  to  a  Mach  disc  reflection 
at  the  appropriate  nozzle  static  pressure  ratio. 

A'~ 

i  \ 

\ 


UMCLAS.STFTF.D 


SECURITY  CL  ASSIFICATIOR  of  Tu,r  PAGEfUTren  Oat*  Erxlrrmd) 


AFWAL-TR-81 -31 61 


FOREWORD 


This  report  is  the  result  of  research  performed  in  the  Computational 
Aerodynamics  Group,  Aerodynamics  and  Airframe  Branch,  Aeromechanics 
Division,  Flight  Dynamics  Laboratory  from  June  1978  to  April  1981.  This 
report  was  prepared  under  Work  Unit  2307N603,  "Computational  Fluid 
Dynamics".  Dr.  Wilbur  L.  Hankey  is  the  Task  Manager. 


The  author  would  like  to  express  his  thanks  to  Dr.  Wilbur  L.  Hankey 
for  his  valuable  advice  and  technical  guidance  in  this  work.  Special 
thanks  are  also  due  to  Dr.  Joseph  S.  Shang  of  the  Computations  I  Aero¬ 
dynamics  Group  for  his  guidance  with  the  numerical  algorithms  portion  of 
the  research.  Being  associated  with  the  members  of  the  Computational 
Aerodynamics  Group  has  been  personally  as  well  as  technically  rewarding. 
Appreciation  is  also  expressed  to  the  members  of  the  AFIT  faculty  who 
served  on  my  doctoral  committee. 

The  material  reported  herein  is  based  on  the  author's  dissertation 
submitted  in  partial  fulfillment  of  the  requirements  for  the  Doctor  of 
Philosophy  degree  at  the  Air  Force  Institute  of  Technology,  Wright- 
Patterson  AFB,  Ohio. 


AFWAL-TR-81 -3161 


TABLE  OF  CONTENTS 

SECTION  PAGE 

I  INTRODUCTION  1 

1 .  Background  1 

2.  Research  Objectives  10 

II  MATHEMATICAL  DESCRIPTION  OF  THE  FLOW  STRUCTURE  12 

1.  Governing  Equations  12 

2.  Boundary  and  Initial  Conditions  19 

III  NUMERICAL  PROCEDURE  23 

1.  Coordinate  System  23 

2.  Solution  Algorithm  27 

3.  Convergence  Criteria  35 

IV  BOUNDARY  AND  INITIAL  CONDITION  IMPLEMENTATION  39 

1 .  The  Upstream  Boundary  39 

2.  The  Upper  Boundary  42 

3.  The  Downstream  Boundary  45 

4.  The  Centerline  46 

5.  The  Nozzle  Walls  48 

6.  Initial  Conditions  50 

V  TURBULENCE  MODELING  52 

1 .  Boundary  Layer  Model  52 

2.  Far  Wake  Model  54 

3.  Near  Wake  Model  59 

VI  NUMERICAL  RESULTS  65 

1 .  Experimental  Data  Base  65 

2.  Computational  Details  66 

3.  Comparison  with  Experimental  Data  73 

VII  CONCLUSIONS  AND  RECOMMENDATIONS  90 

APPENDIX  A  NOZZLE  WALL  TEMPERATURE  CALCULATION  93 

APPENDIX  B  ADAPTIVE  FINITE  DIFFERENCE  MESH  96 

APPENDIX  C  TWO-DIMENSIONAL  FLAT  PLATE  FAR  WAKE  SOLUTION  99 

APPENDIX  D  TWO-DIMENSIONAL  WEDGE-FLAT  PLATE  NEAR  WAKE 

SOLUTION  103 

APPENDIX  E  AXISYMMETRIC  NOZZLE  SOLUTION  SIMULATING 

INTERNAL  BOUNDARY  LAYER  SEPARATION  110 

APPENDIX  F  INVESTIGATION  OF  NUMERICAL  ERROR  114 

REFERENCES  120 


v 


AFWAL-TR-81 -31 61 


LIST  OF  ILLUSTRATIONS 

FIGURE  PAGE 

1  Shock  Structure  for  a  Typical  Overexpanded  Axisymmetric 

Nozzle  4 

2  Schematic  of  a  Typical  Mach  Reflection  5 

3  Hodograph  Diagram  Showing  Shock  Polar  Intersections 

for  Transition  from  Regular  to  Mach  Reflection  7 

4  Viscous  Effects  on  the  Mach  Disc  Structure  9 

5  Schematic  of  the  Axisymmetric  Jet  Model  Used  as  a 

Basis  for  the  Computational  Solutions  11 

6  Physic-'  Domain  for  the  Computational  Solutions  19 

7  Finite  Difference  Mesh  in  Physical  Space  24 

8  Exponentially  Stretched  Mesh  Schematic  24 

9  Adaptive  Finite  Difference  Mesh,  P./P  =  0.150  26 

J  00 

10  Computational  Mesh  in  the  Transformed  Plane  26 

11  Graphical  Representation  of  the  Numerical  Sweep 

Operators  30 

12  Base  Pressure  Tap  Location  in  the  Flowfield  37 

13  Typical  Base  Pressure  Convergence,  P./P^  =  0.527  38 

14  Flowfield  Schematic  for  an  Axisymmetric  Nozzle  in  a 

Supersonic  External  Stream  40 

15  Computed  Static  Pressure  Variation  Along  the  Ogive 
Body  Using  a  Parabolized  Navier-Stokes  Solver 

(Reference  23)  41 

16  Upper  Boundary  Condition  Schematic  43 

17  Finite  Difference  Mesh  Near  the  Nozzle  Walls  49 

18  Eddy  Viscosity  Model  Domains  53 

19  Eddy  Viscosity  Distribution  Across  a  Boundary  Layer  65 

20  Typical  Vorticity  Profile  Used  to  Compute  Mixing 

Layer  Thickness  57 


vi 


AFWAL-TR 

-81-3161 

LIST  OF  ILLUSTRATIONS  (Cont'd) 

FIGURE 

PAGE  ' 

21 

Maximum  Velocity  Defect  vs  Distance  Behind  the  Trailing 

Edge  of  a  Two-Dimensional  Flat  Plate,  =  1.60 

58 

22 

Viscous  Layer  Structure  in  the  Near  Wake 

60 

23 

Computed  Nozzle  Base  Pressure  vs  the  Position  of  the 

Mixing  Region  Midpoint  Xm 

62 

24 

Computed  Mixing  Layer  Thickness  Used  in  the  Wake 

Turbulence  Model  for  the  Wedge-Flat  Plate  Case 

63  j 

25 

Adaptive  Finite  Difference  Mesh,  P./P  =  0.527 

3  00 

68  | 

26 

Error  in  Computed  Velocity  Gradient  vs  Grid  Courseness 
for  a  Turbulent  Boundary  Layer 

70 

27 

Axisymmetric  Nozzle  Solution,  P./P^  *  0.150 

74 

28 

Axisymmetric  Nozzle  Solution,  Pj/P^  =  0.251 

75 

|  29 

Axisymmetric  Nozzle  Solution,  PJPm  =  0.527 

~  76 

30 

Axisymmetric  Nozzle  Solution,  P./P^  =  1-03 

77 

« 

31 

Axisymmetric  Nozzle  Solution,  Pj/P^  =  1.59 

78 

32 

Computed  Mach  Number  Contours  in  the  Region  Near  the 

Shock  Reflection  at  the  Centerline 

79 

33 

Computed  Velocity  Profiles,  P./Pm  =  0. i 50 

81 

34 

Computed  Velocity  Fields  in  the  Near  Wake  Region  of 
the  Nozzle  Base  Annulus 

82 

35 

Position  of  the  Dividing  Streamline  in  the  Computational 
Nozzle  Solutions 

83 

36 

Axial  Variation  in  Computed  Centerline  Mach  Number, 

P./P  =  0.150 

J  00 

84 

37 

Shock  Reflection  Lengths  Along  the  Nozzle  Centerline 
vs  Nozzle  Pressure  Ratio 

85 

38 

Base  Pressure  Coefficient  vs  Nozzle  Pressure  Ratio 

87 

39 

Schlieren  Photographs  (Reference  16)  Showing  the  Eventual 
Deterioration  of  the  Mach  Disc  with  Decreasing  Nozzle 
Pressure  Ratio 

i  ; 

89  . 

* 

vi  i 

t 

_ j 

AFWAL-TR-81 -3161 


LIST  OF  ILLUSTRATIONS  (Concluded) 


FIGURE 

PAGE 

40 

Heat  Flux  Balance  Used  to  Determine  the  Nozzle  Wall 
Temperature 

94 

41 

Adaptive  Mesh  Schematic 

98 

42 

Velocity  Profile  at  the  Trailing  Edge  of  the  Two- 
Dimensional  Flat  Plate 

100 

43 

Computational  Mesh  Used  in  the  Flat  Plate  Solution 

102 

44 

Computed  Flat  Plate  Velocity  Profiles 

102 

45 

Computational  Boundary  Conditions  for  the  Two- 
Dimensional  Wedge-Flat  Plate 

104 

46 

Computational  Mesh  Used  for  the  Two-Dimensional 

Wedge-Flat  Plate  Solution 

106 

47 

Two-Dimensional  Wedge-Flat  Plate  with  the  Computed 

Mach  Number  Contours  Shown 

107 

48 

Computed  Velocity  Profiles  in  the  Near  Wake  of  the 
Two-Dimensional  Wedge-Flat  Plate 

108 

49 

Static  Pressure  Along  the  Line  of  Symmetry  in  the 

Near  Wake  of  the  Two-Dimensional  Wedge-Flat  Plate 

108 

50 

Pitot  Pressure  Profiles  in  the  Near  Wake  of  the  Two- 
Dimensional  Wedge-Flat  Plate  (Symbol s-Experimental  Data 
(Reference  34),  Solid  Lines  -  Computational  Solution) 

109 

51 

Computed  Velocity  Profiles  Near  the  Nozzle  Annulus  for 
the  Separated  Flow  Simulation 

112 

52 

Computed  Mach  Number  Contours  for  the  Separated  Flow 
Simulation 

113 

53 

Comparison  between  the  Computational  Solution  and  the 
Error  Present  in  the  Continuity  Equation,  P./P  =  0.150 

J  CO 

116 

54 

Extension  of  the  Downstream  Boundary  Showing  Computed 

Mach  Number  Contours,  P./P  =  0.251 

J  “ 

118 

55 

Extension  of  the  Upper  Boundary  Showing  Computed  Mach 

Number  Contours,  P./P  =  0.251 

J  05 

119 

AFWAL-TR-81-3161 


LIST  OF  TABLES 


TABLE 

PAGE 

1 

Computational  Grid  Parameters 

67 

2 

Computational  Jet  Flow  Parameters 

71 

3 

Comparison  Between  a  1-D  Analysis  and  the 
Solution  Across  the  Mach  Disc  for  Pj/P^  = 

Computational 

0.150 

80 

4 

Comparison  of  Shock  Reflection  Lengths 

86 

5 

Comparison  of  Mach  Disc  Radii 

86 

6 

Comparison  of  Base  Pressure  Coefficients 

87 

AFWAL-TR-81 -31 61 


e 

F 

FC 

f 

G 

GC 

H 

H 


hi 


IL 

IW 

JL 

JWI 

JWO 

j 


LIST  OF  SYMBOLS 

Constant  used  to  obtain  exponentially  stretched  mesh. 
Adaptive  mesh  constant. 

2 

Local  skin  friction  coefficient,  2x  /p  u  . 

Specific  heat  at  constant  pressure. 

Specific  heat  at  constant  volume. 

Speed  of  sound,  /yRT. 

Error  vector. 

Energy,  CyT  +  (u2  +  v2)/2. 

Flux  vector.  Equation  28. 

Damping  vector.  Equation  30. 

Primitive  flow  variable  p,  u,  v,  or  T. 

Flux  vector.  Equation  28. 

Damping  vector.  Equation  32. 

Flux  vector.  Equation  28 

Total  Enthalpy,  CpT  +  (u2  +  v2)/2. 

Convective  heat  transfer  coefficient,  Appendix  A. 

Thickness  of  the  two-dimensional  wedge-flat  plate. 

Index  for  grid  points  in  the  axial  direction. 

Total  number  of  grid  points  in  the  axial  direction. 

Number  of  grid  points  axially  along  the  nozzle  wall. 

Total  number  of  grid  points  in  the  radial  direction. 

Radial  index  of  the  inner  nozzle  wall. 

Radial  index  of  the  outer  nozzle  wall. 

Index  for  grid  points  in  the  radial  direction. 

Exponent  parameter  equal  to  0  or  1  for  two-dimensional 
or  axisymmetric  flow,  respectively. 


x 


AFWAl-TR-81 -3161 


K 


K 


o 


k 

L 


L 


b 


L 


m 


L 


l 


M 


PB 

Pr 

-> 

q 

V  qr 

R 

Re 

r 

r 

gr 

rm 

S 


T 

t 


LIST  OF  SYMBOLS  (Cont'd) 

Ratio  of  radial  grid  point  heights,  r ( i ,3)/r(i ,2) . 

Initial  velocity  profile  constant  in  the  wake  region. 
Equation  111. 

Thermal  conductivity. 

Length  scale. 

Experimental  ogive  body  length 

Length  scale  used  to  generate  stretched  mesh. 

MacCormack  difference  operators  in  the  n  and  £  directions. 
Prandtl  turbulent  mixing  length. 

Mach  number,  (u^  +  v^)^/c. 

Static  pressure. 

Base  pressure. 

Pitot  pressure. 

Base  pressure  coefficient,  (Pg  -  P^J/q^. 

Prandtl  number,  uCp/k. 

Heat  flux  vector. 

Heat  flux  in  the  axial  and  radial  directions. 

Gas  constant  for  air. 

2  2  1/2 

Reynolds  number,  p(u  +  v  )  '  L/y. 

Spatial  coordinate  normal  to  the  nozzle  centerline. 

Length  of  the  computational  grid  in  the  radia1  direction. 
Mach  disc  radius. 

Fluid  stress  tensor  that  includes  pressure  and  viscous 
forces. 

Absolute  temperature. 

Time. 


xi 


cl 


AFWAL-TR-81 -3161 


LIST  OF  SYMBOLS  (Cont'd) 


t 


ch 


U 


u 

uch 

+ 

u 

V 

v 

X 


X 


r 


x 


gr 


y 


yB 

+ 

y 

GREEK  SYMBOLS: 

V  as 

Y 

6 

6 

6* 

6ij 

E 


n 

A 

0 

0 


Characteristic  time.  Equation  84. 

Conservative  flow  variable  vector.  Equation  28. 

Velocity  vector. 

Velocity  component  along  the  x  axis. 

Characteristic  velocity.  Equation  84. 

1 12 

Nondimensional  velocity  component,  u/(t  /p) 

Volume. 

Velocity  component  along  the  r  axis. 

Spatial  coordinate  parallel  to  the  nozzle  centerline. 
Mach  disc  reflection  length. 

Length  of  the  computational  grid  in  the  axial  direction 

Spatial  coordinate  normal  to  the  x  axis  in  a  two- 
dimensional  flow. 

Radial  thickness  of  the  nozzle  base  annulus. 

1  /2 

Nondimensional  height,  y (t  /p )  /v. 


Damping  coefficients  in  the  n  and  c  directions. 

Ratio  of  specific  heats,  Cp/ Cy . 

Flow  deflection  angle.  Section  I. 

Momentum  boundary  layer  thickness. 

Momentum  boundary  layer  displacement  thickness. 
Kronecker  delta,  equal  to  0  if  i/j  or  equal  to  1  if  i=j 
Turbulent  eddy  viscosity  coefficient. 

Transformed  coordinate  normal  to  the  nozzle  centerline. 
Designates  a  finite  difference  when  used  as  a  prefix. 
Local  flow  angle,  arctan(v/u) 

Boundary  layer  momentum  thickness,  Appendix  C. 


xi  i 


AFWAL-TR-81-3161 


A 

XI 

M 

^M 

v 

? 

P 

a 

Txr 

0) 

SUBSCRIPTS: 

aw 

e 

i » J 
i  >  j 
j 

n 

s 

t 

tr 


w 

0 


00 


LIST  OF  SYMBOLS  (Cont’d) 

Viscosity  diffusion  coefficient.  Equation  37. 

Left  running  characteristic  line. 

Absolute  viscosity  coefficient. 

Mach  angle  of  supersonic  flow.  Equation  49. 

Kinematic  viscosity,  p/p. 

Transformed  coordinate  parallel  to  the  nozzle  centerline. 
Fluid  density. 

Normal  stress  on  an  element  of  fluid.  Equations  29-31. 
Shear  stress.  Equation  32. 

Vorticity. 

Adiabatic  wall. 

Evaluated  at  the  edge  of  the  boundary  layer. 

Grid  point  indices. 

Indicial  notation  Section  II. 

Evaluated  in  the  jet  flow. 

Normal  to  the  nozzle  wall. 

Tangent  to  the  nozzle  wall. 

Turbulent  flow. 

Transition  from  regular  shock  reflection  to  Mach  disc 
reflection. 

Evaluated  in  the  wake. 

Stagnation  value. 

Evaluated  in  the  external  freestream. 


xiii 


AFWAL-TR-81 -31 61 


SUPERSCRIPTS: 

n 

n+1/2 

n+1 

o 


LIST  OF  SYMBOLS  (Concluded) 

Evaluated  at  known  time,  nAt. 

Evaluated  at  intermediate  predictor  time  level, 
(n+l/2)At. 

Evaluated  at  new  corrector  time  level,  (n+1 ) A t . 
Evaluated  at  time  level  t=0. 


OTHER  NOTATION: 

(=)  Oenotes  a  two-dimensional  matrix. 

(")  Oenotes  time  averaged  values.  Section  II. 

(  )■  Oenotes  unsteady  values  due  to  turbulence.  Section  II. 


xi  v 


AFWAL-TR-81 -31 61 


SUMMARY 

The  use  of  computational  analysis  in  the  design  of  propulsive  nozzle 
installations  has  recently  expanded  as  advanced  digital  computers  have 
been  developed  which  result  in  lowering  computational  costs  versus 
actual  wind  tunnel  test  costs.  Although  a  range  of  numerical  techniques 
has  been  applied  in  this  area,  only  those  utilizing  the  full  Navier- 
Stokes  equations  across  the  flow  domain  have  successfully  simulated  the 
viscous  phenomena  associated  with  aft-end  flowfields.  Navier-Stokes 
methods  are  particularly  useful  for  predicting  off-design  nozzle 
characteristics  where  the  overexpanded  or  underexpanded  flowfield  is 
more  complex  and  where  viscous  regions  are  more  prevalent  than  at  on- 
design  conditions.  One  feature  typical  of  these  off-design  conditions 
is  the  appearance  of  a  strong  normal  shock  wave  referred  to  as  a  Mach 
disc.  Viscous  nozzle  flowfields  containing  this  phenomenon  have  not  been 
adequately  simulated  in  the  past.  This  research  details  the  development 
of  a  numerical  Navier-Stokes  method  capable  of  accurately  predicting 
nozzle  flowfields  which  contain  both  highly  viscous  regions  and  complex 
shock  structures  typified  by  the  Mach  disc  formation. 

Numerical  solutions  to  the  Navier-Stokes  equations  are  obtained  for 
a  domain  containing  an  axisymmetric  nozzle  in  a  supersonic  external 
stream  (M^  =  1.94,  Mjet  =  3.0,  Re^  =  2.2x10^).  Five  nozzle  pressure 

ratio  conditions  ranging  from  a  highly  overexpanded  case  ( P  - / P  =  0.15) 

J  00 

which  exhibits  a  Mach  disc  shock  formation,  to  a  slightly  underexpanded 
case  (P./P^  =  1.59)  are  examined  and  solved  numerically.  The  weak  con- 
servative  form  of  the  two-dimensional  (axisymmetric),  time  dependent 
Navier-Stokes  equations  is  solved  using  MacCormack's  explicit  finite 
difference  method.  This  algorithm  is  an  efficient  Lax-Wendroff  type 
differencing  scheme  of  second  order  accuracy  which  utilizes  time¬ 
splitting  and  two-step  predictor-corrector  techniques.  An  adaptive 
grid  scheme  is  utilized  in  the  wake  of  the  nozzle  base  annulus  that 
allows  the  fine  mesh  region  of  the  computational  grid  to  remain  in  the 
mixing  layer  containing  high  flow  gradients  as  each  solution  progresses 
towards  convergence.  Appropriate  numerical  boundary  conditions  are 
applied  that  allow  the  computational  domain  to  be  restricted  to  a 
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compact  region  surrounding  the  nozzle.  Locally  dependent  eddy  viscosity 
modelling  is  applied  in  the  form  of  a  Cebeci-Smith  two-layer  model  in  the 
boundary  layer  regions  on  the  nozzle  walls,  and  a  form  of  the  Prandtl 
mixing  length  model  in  the  nozzle  wake  region. 

The  numerical  solutions  successfully  reproduced  all  of  the  essential 
nozzle  flow  features  including  boundary  layers,  corner  expansions, 
recompression  shocks,  the  separated  recirculation  region  along  the 
nozzle  base  wall,  and  the  evolution  of  the  near  wake  to  a  far  wake  type 
of  flow.  Correct  transition  from  regularly  reflected  shock  waves  at  the 
line  of  symmetry  in  the  jet  core  to  the  strong  Mach  disc  shock  reflection 
was  numerically  achieved,  as  was  the  simulation  of  the  subsonic  embedded 
region  immediately  behind  the  Mach  disc  shock  structure.  Numerically 
obtained  nozzle  base  pressure  coefficients  were  within  seven  percent  of 
the  experimentally  determined  values  for  all  cases  where  the  attached 
flow  assumption  in  the  divergent  portion  of  the  nozzle  was  valid. 

Present  solution  times  for  2,500  point  grids  are  on  the  order  of 
two  to  four  hours  when  run  on  a  Cyber  750  computer.  A  fully  vectorized 
version  of  the  present  computer  code  can  be  expected  to  converge  within 
five  minutes  on  a  CRAY-1  computer  for  similar  grids,  allowing  the  compu¬ 
tation  of  more  complex  nozzle  geometries  and  better  resolution  in  the 
boundary  layers  through  the  use  of  a  finer  mesh  in  future  efforts. 
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SECTION  I 
INTRODUCTION 

1 .  BACKGROUND 

The  increased  importance  of  the  aft-end  drag  problem  associated  with 
nozzle  installations  in  current  and  future  high  performance  aircraft  has 
led  to  extensive  and  very  costly  experimental  nozzle  test  programs. 

Any  technique  which  can  reduce  this  requirement  for  wind  tunnel  testing 
in  the  design  of  nozzle  installations  will  result  in  a  significant 
savings  to  the  technical  community  of  both  time  and  resources. 

Computational  aerodynamics  shows  great  promise  as  a  field  which 
can  have  a  favorable  impact  on  this  requirement  for  nozzle  design 
information.  Current  computational  techniques  in  this  area  utilize 
advanced  digital  computers  to  simulate  the  flowfield  surrounding  the 
nozzle  at  projected  flight  conditions.  It  has  been  shown  that  boundary 
layer  and  shear  layer  growth,  areas  of  separated  flow,  shock  wave 
formation  and  interactions,  and  jet  plume  blockage  and  entrainment 
characteristic  of  nozzle  flows  can  be  analyzed  using  computational 
techniques.  Unlike  experimental  testing,  computational  analysis  is 
not  necessarily  restricted  by  wind  tunnel  Reynolds  number  or  nozzle 
exhaust  temperature  limitations.  Flowfields  analyzed  computationally 
can  also  eliminate  the  undesirable  effects  of  support  stings  and  test 
section  walls  that  occur  during  experimental  testing.  As  more  advanced 
computers  are  developed,  the  cost  of  numerically  analyzing  a  given 
nozzle  configuration  is  decreasing.  Since  the  cost  of  wind  tunnel 
testing  is  steadily  increasing,  computational  analysis  is  being 
utilized  more  extensively. 

Several  of  the  first  computational  solutions  to  include  viscous 
effects  inherent  to  aft  end  or  nozzle  flowfields  consisted  of  patching 
techniques  that  divided  the  field  into  predominantly  inviscid  and 
viscous  regions.  Grossman  and  Melnik  (Reference  1),  and  Cosner  and 
Bower  (Reference  2)  obtained  transonic  boattail  nozzle  solutions  using 
iterative  techniques  that  divided  the  flowfield  into  an  inviscid  free- 
stream,  an  inviscid  jet,  and  a  viscous  boundary  layer  and  mixing  layer 
region.  The  freestream  solution  procedure  assumed  irrotational  potential 
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flow  that  could  be  solved  by  a  relaxation  algorithm  applied  to  the 
potential  flow  equations.  The  rotational  inviscid  supersonic  jet  was 
solved  using  a  hyperbolic  marching  technique.  Imbedded  shocks  in  the 
jet  were  explicitly  fitted  to  satisfy  the  Rankine-Hugoniot  equations. 

The  viscous  mixing  region  was  assumed  to  be  isobaric  and  was  solved  using 
integral  techniques.  Each  region  was  solved  separately  and  patched 
together  iteratively  using  pressure  and  flow  direction  conditions  at  the 
common  boundaries.  Separation  regions  could  not  be  accounted  for,  so 
equivalent  fitted  body  blending  was  used  to  obtain  reasonable  flow 
solutions. 

Pergament,  Dash,  and  Wilmoth  (Reference  3)  introduced  a  displacement 
thickness  correction  to  the  inviscid  plume  boundary  to  account  for  the 
effects  of  jet  entrainment  on  the  inviscid  external  flow  calculation. 
Their  analysis  also  included  the  effects  of  species  mixing  and  pressure 
gradients  in  the  mixing  region,  but  still  could  only  account  for 
separation  regions  by  body  blending  techniques.  Yeager  (Reference  4) 
attempted  to  include  a  fourth  separation  region  involving  reci rculating 
flow  that  was  defined  by  a  dividing  streamline  which  connected  separation 
and  reattachment  points.  The  extent  of  this  region  was  determined  using 
local  control  volume  analyses,  and  it  was  found  that  reasonable 
reattachment  points  could  only  be  predicted  through  the  application  of 
empirical  corrections  during  the  solution  procedure.  Although  these 
iteratively  patched  solutions  gave  reasonable  results  for  specific  data 
sets,  the  required  amount  of  empirical  matching  and  explicit  fitting 
limits  the  use  of  this  type  of  computational  method  as  a  predictive 
technique. 

A  more  adaptable  method  of  simulating  the  viscous-invi scid  inter¬ 
actions  that  occur  in  typical  nozzle  flowfields  involves  solving  the 
time  dependent,  compressible  Navier-Stokes  equations  uniformly  over  the 
entire  nozzle  flowfield.  This  approach  has  a  direct  advantage  over  the 
previously  discussed  iteratively  patched  methods  where  an  accurate 
viscous-inviscid  matching  procedure  is  required  in  order  to  obtain 
reasonable  results.  In  the  direct  approach,  the  predominantly  inviscid 
and  viscous  flow  regions  are  computed  simultaneously  with  no  matching 
required.  Holst  (Reference  5)  used  this  approach  to  solve  for  supersonic 
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flow  over  axisymmetric  boattail  nozzles  with  plume  simulators.  Although 
a  plume  simulator  does  not  model  the  entrainment  and  blockage  of  a  jet 
plume,  its  flowfield  does  contain  phenomena  characteristic  to  nozzles 
such  as  turbulent  boundary  layers,  recompression  shock  waves,  and  separated 
recirculating  regions  of  flow.  Holst's  solutions  were  obtained  using 
MacCormack's  explicit  finite  difference  algorithm,  a  stretched  mesh 
aligned  with  the  solid  body  through  an  analytic  transformation  and  a 
two-layer  eddy  viscosity  model  to  account  for  the  Reynold's  stresses 
that  included  a  relaxation  formula  to  model  the  separated  flow  region. 
Pressure  distributions,  skin  friction  coefficients  and  areas  of  separated 
flow  were  in  good  agreement  with  experimental  data,  particularly  in  the 
cases  where  a  fine  mesh  was  utilized.  Mikhail  (Reference  6)  recently 
computed  solutions  for  viscous  supersonic  flow  around  an  axisymmetric 
boattail  nozzle  with  a  jet  exhaust  flow.  MacCormack's  explicit  method 
was  again  used  as  the  numerical  algorithm,  together  with  a  surface 
oriented  mesh  system  obtained  through  a  numerical  mapping  procedure. 
Reynold's  stresses  were  also  accounted  for  through  the  application  of 
algebraic  eddy  viscosity  models.  Reasonable  agreement  with  experimental 
surface  pressure  data  on  the  boattail  was  obtained. 

Navier-Stokes  solutions  are  especially  useful  for  predicting  the 
off-design  nozzle  characteristics  where  the  flowfield  is  in  either  a 
significantly  overexpanded  or  underexpanded  state.  At  these  conditions 
the  flow  structure  is  usually  more  complex  with  viscous  regions  becoming 
more  prevalent  than  at  on-design  conditions.  One  feature  typical  of 
these  off-design  conditions  is  the  establishment  of  a  triple-point  in 
the  jet  flow,  and  the  appearance  of  a  strong  normal  shock  wave  referred 
to  as  a  Mach  disc  in  axisymmetric  flow  or  a  Riemann  wave  in  two-dimensional 
flow  (Figure  1).  This  strong  shock  formation  occurs  when  the  deflection 
angle  of  the  jet.  flow  is  large  enough  so  that  the  resulting  shock  wave  is 
too  strong  for  a  regular  reflection  at  the  centerline  to  exist.  Near  the 
centerline  the  Mach  disc  must  be  normal,  since  this  is  the  only  way  a 
shock  can  occur  without  any  change  in  flow  direction.  As  shown  in 
Figure  2,  both  the  Mach  disc  and  the  reflected  shock  are  curved  near 
the  triple  point  (Reference  7).  A  slip  line  emanates  from  the  triple 
point,  and  the  flow  downstream  of  the  Mach  disc  and  reflected  shock  is 
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MACH  REFLECTION 


Figure  1.  Shock  Structure  for  a  Typical  Overexpanded  Axisymmetri 
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REFLECTED 


Figure  2.  Schematic  of  a  Typical  Mach  Reflection 
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rotational  in  nature  due  to  the  curvature  of  the  shocks.  As  discussed  by 
Henderson  and  Lozzi  (Reference  8),  this  region  downstream  of  the  shocks 
may  be  either  totally  subsonic  or  contain  both  supersonic  and  subsonic 
regions.  If  the  incident  Mach  number  M-]  is  greater  than  2.40,  region  3 
will  be  supersonic,  while  region  4  remains  subsonic.  Since  the  jet  Mach 
numbers  relevant  to  this  investigation  are  greater  than  2.40,  the  latter 
case  in  which  a  subsonic  core  only  exists  behind  the  Mach  disc  will  be 
examined. 

The  transition  from  regular  to  Mach  reflections  can  also  be 
examined  using  a  hodograph  diagram  shown  in  Figure  3.  For  deflection 
angles  less  than  the  transition  angle  the  flow  can  be  brought  to 

the  required  zero  deflection  at  state  3  through  a  weak  regular  wave 
reflection.  For  deflection  angles  greater  than  the  transition  angle,  the 
flow  in  region  3  cannot  achieve  a  zero  deflection  state,  and  lies  on  the 
strong  shock  portion  of  the  initial  shock  polar.  The  flow  near  the 
centerline  passes  through  a  strong  normal  shock  to  condition  4'  with  no 
deflection  occurring.  The  flow  state  at  the  curved  Mach  disc  then 
exists  along  the  strong  shock  portion  of  the  incident  shock  polar  from  a 
zero  deflection  state  4'  near  the  centerline  to  a  deflected  state  4  near 
the  triple  point  with  a  pressure  and  flow  angle  equal  to  that  in  region  3, 
but  with  different  velocity  and  entropy  values  that  generate  the  slip 
line. 

The  mixed  supersonic-subsonic  flow  region  surrounding  the  Mach  disc 
greatly  complicates  the  analysis  of  nozzle  flowfields  in  which  this  shock 
structure  is  present.  Flowfields  containing  this  phenomena  have  not  been 
adequately  simulated  in  the  past  using  viscous  techniques.  A  variety  of 
techniques  for  locating  the  triple  point  and  the  resulting  normal  shock 
have  been  presented  that  utilize  an  iterative  combination  of  the  method 
of  characteristics  and  schemes  involving  approximate  analyses  such  as 
pressure  requirements  downstream  of  the  strong  shock  (References  9  and  10) 
or  one-dimensional  flow  calculations  downstream  through  a  throat  region 
in  the  flow  (References  11  and  12).  Although  each  of  these  methods  give 
reasonable  results  for  determining  the  triple  point  location  and  size  of 
the  resulting  Mach  disc,  each  is  only  valid  for  a  limited  range  of  nozzle 
pressure  ratios  and  jet  Mach  numbers.  In  addition,  none  of  these  techniques 
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give  a  solution  for  the  subsonic  core  region  downstream  of  the  normal 
shock.  At  least  two  time-dependent  inviscid  techniques  have  been  used 
to  overcome  the  deficiencies  of  the  semi -empirical  methods  previously 
mentioned.  Jofre  (Reference  13)  performed  a  finite  difference  technique 
for  an  underexpanded  jet  with  a  Mach  disc  solution.  A  method  of  char¬ 
acteristics  solution  was  used  in  the  plume  expansion  region  near  the 
nozzle  exit.  This  gave  an  upstream  flow  profile  used  in  the  time 
dependent  solution  further  downstream.  Sinha,  Zakkay,  and  Erdos 
(Reference  14)  analyzed  a  two-dimensional  underexpanded  jet  containing  a 
strong  normal  shock  using  a  finite  difference  technique  over  the  entire 
flowfield  of  interest.  Both  of  these  investigations  used  versions  of 
Lax-Wendroff  numerical  algorithms  and  simple  square  grids.  These 
solutions  were  much  more  adaptable  than  the  previous  semi-empirical 
techniques  since  the  flow  tends  to  adjust  to  its  local  environment  so 
that  the  proper  shock  structure  is  automatically  obtained  as  the  solution 
develops.  However,  these  solutions  represent  only  a  first  approximation 
of  the  correct  viscous  solution  (Reference  15),  particularly  in  the  region 
downstream  of  the  normal  shock  as  shown  in  Figure  4.  Flow  properties  in 
this  region  were  found  to  be  heavily  dependent  on  the  level  of  damping 
used  in  the  solution  algorithm.  For  example,  Jofre  (Reference  13)  found 
that  an  unrealistic  region  of  reverse  flow  was  generated  immediately 
behind  the  normal  shock  unless  heavy  damping  was  applied  in  the  solution 
procedure. 

All  of  these  inviscid  solutions  involved  jets  exhausting  into  a 
quiescent  atmosphere.  Solutions  were  not  attempted  for  the  more  complex 
case  of  a  nozzle  in  which  the  external  flow  stream  interacts  with  the 
jet.  A  full  Navier-Stokes  solution  which  accounts  for  the  viscous 
effects  present  in  the  flowfield  at  these  off-design  conditions  is 
necessary  in  order  to  adequately  simulate  both  the  strong-shock  structure 
with  its  resulting  imbedded  subsonic  flow  region  as  well  as  the  inter¬ 
action  of  the  jet  plume  with  the  external  flowfield. 
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2.  RESEARCH  OBJECTIVES 

The  primary  objective  of  this  research  is  the  development  of  a 
numerical  Navier-Stokes  method  capable  of  accurately  predicting  nozzle 
flowfields  which  contain  both  highly  viscous  regions  and  complex  shock 
structure  typified  by  the  Mach  disc  shock  formation.  Overexpanded 
axi symmetric  nozzles  will  primarily  be  simulated,  since  they  meet  the 
previous  criteria  while  possessing  fairly  compact  flow  domains  which 
contain  the  flow  phenomena  of  interest.  The  experimental  data  of  Bromm 
and  O'Donnell  (Reference  16)  has  been  chosen  as  a  basis  for  comparison 
in  this  research  effort.  Data  in  this  reference  is  given  for  an 
axisymmetric  Mach  three  isentropic  nozzle  embedded  in  a  turbulent 
Mach  1.94  external  flowfield  as  shown  in  Figure  5.  Nozzle  pressure  ratios 
ranging  from  a  slightly  underexpanded  condition  to  a  highly  overexpanded 
condition  which  exhibits  the  Mach  disc  structure  were  obtained  experi¬ 
mentally.  This  particular  nozzle  possesses  a  relatively  thick  base 
annulus  which  generates  a  strong  vi scous-invi scid  interaction  in  the 
near  wake  region  of  the  nozzle.  These  interactions  affect  the  develop¬ 
ment  of  the  primarily  inviscid  shock  structure,  and  can  only  be  analyzed 
properly  using  a  full  Navier-Stokes  methodology 

Although  Mikhail  achieved  full  Navier-Stokes  solutions  for  an 
axisymmetric  boattail  nozzle,  he  was  not  able  to  generate  an  accurate 
solution  for  the  condition  at  which  a  Mach  disc  shock  structure  was 
shown  to  exist  experimentally  (Reference  6).  Possible  causes  of  this 
inability  to  generate  the  strong  shock  structure  include  boundary 
condition  formulation,  mesh  spacing  and  turbulence  modeling.  These 
three  areas  will  be  concentrated  on  in  the  present  investigation  in 
order  to  achieve  the  desired  goal  of  an  accurate  predictive  technique  for 
both  on-design  and  off-design  nozzle  performance. 
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SECTION  II 

MATHEMATICAL  DESCRIPTION  OF  THE  FLOW  STRUCTURE 
1.  GOVERNING  EQUATIONS 

The  governing  equations  for  flows  containing  the  shock  and  viscous 
phenomena  of  interest  are  the  conservation  equations  for  mass,  momentum, 
and  energy  known  as  the  Navier-Stokes  equations.  The  gases  involved  are 
assumed  to  be  single  component,  have  constant  specific  heats,  and  obey 
the  perfect  gas  equation  of  state: 

P  -  pRT  (1 ) 


In  computational  fluid  dynamics,  the  Eulerian  method  is  usually 
applied  to  the  problem  of  interest.  This  method  involves  a  fixed  control 
volume  that  is  specified  relative  to  a  given  coordinate  system.  Properties 
of  the  fluid  are  then  specified  as  functions  of  both  time  an"1  space.  The 
conservation  equations  are  approached  using  this  methodology. 

a.  Conservation  of  Mass 


For  a  given  system  in  which  matter  is  neither  created  or  destroyed, 
the  law  of  mass  conservation  can  be  written  as 


+  V  •  pu)  dV 


0 


(2) 


where  V  is  an  arbitrary  volume  fixed  in  space, 
b.  Conservation  of  Momentum 

For  a  given  system,  the  law  of  momentum  conservation  states  that 
the  rate  of  change  of  momentum  is  equal  to  the  sum  of  the  external  forces 
acting  on  the  control  volume.  If  body  forces  are  neglected,  this  law 
can  be  written  as: 

+  v  •  («>  I]  d* 
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The  variable  S  denotes  a  stress  tensor  involving  pressure  and  viscous 
forces  which  acts  on  the  fluid. 

c.  Conservation  of  Energy 

The  law  of  conservation  of  energy  states  that  for  a  given  system 
which  does  not  contain  any  internal  heat  sources,  the  rate  of  change 
of  the  total  energy  of  the  system  is  equal  to  the  heat  added  into  the 
system  plus  the  work  done  on  the  system  by  viscous  and  pressure  forces. 
This  can  be  stated  as 


///. 


1  3t 


(pe)  u]  dV 


III 


(V-u*S  -  V *q)  dV 


(4) 


Since  these  conservation  equations  are  valid  for  any  arbitrary 
volume  V;  when  the  integrands  are  continuous,  these  equations  imply  that: 

3p/3t  +  v  «  pu  •  0  (5) 

3(pu)/3t  +  V  •  (pu)  u  -  V  •  S  »  0  (6) 

3(pe)/3t  +  V  •  (pe)  u  +  (V*q  -  V*u’S)  ■  0  (7) 

It  should  be  noted  that  these  equations  are  written  in  conservative 
form  where,  for  the  two-dimensional  and  axisymmetric  flows  of  interest, 
the  applicable  dependent  variables  are  p,  pu,  pv,  and  pe.  As  shown  by 
Roache  (Reference  17),  this  conservative  form  allows  the  finite  difference 
equations  to  preserve  the  Gauss  divergence  property  of  the  continuum 
equations.  This  form  allows  a  oalance  between  the  flux  quantities  and 
accumulation  rates  for  a  small  control  volume.  Roache  also  states  that 
the  Rankine-Hugoniot  shock  relations  were  derived  using  the  conservative 
form.  Thus,  shock  jump  conditions  are  automatically  satisfied  since  the 
conservative  variables  are  continuous  across  the  shock  and  need  no 
special  treatment  because  of  discontinuities.  This  approach  is  known  as 
shock  capturing  or  shock  smearing.  The  conservation  form  of  the  equations 
then  allows  the  finite  difference  formulation  to  satisfy  the  physical  laws 
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on  a  macroscopic  scale,  not  merely  in  some  academic  limit  as  Ax,  Ay, 
and  At  approach  zero. 

Since  the  flowfields  of  interest  are  turbulent,  the  solution  of  the 
conservation  equations  must  take  into  account  the  effects  of  the  random 
fluctuations  of  the  dependent  variables  inherent  to  turbulent  flows.  In 
accounting  for  these  effects,  Cartesian  tensor  notation  will  be  applied. 
The  usual  conventions  of  a  repeated  subscript  indicating  summation  over 
the  entire  range  of  indices  and  a  comma  representing  partial  differentia¬ 
tion  will  be  used  to  make  the  equations  compact.  Cartesian  tensors  are 
used  to  allow  working  directly  with  the  physical  components,  while  still 
being  applicable  to  the  2-D  and  axi symmetric  systems  of  interest.  The 
conservation  Equations  5  through  7  can  then  be  written  as: 


p,t  +  “  0 

(pui>’t  +  (pUiUj  +  6ij  P 
(pe) , t  +  (Peuj  +  “  ui 


'  V-j 

V’j  ‘ 


0 


0 


(8) 

(9) 

(10) 


where  the  stress  tensor  has  been  expanded  in  the  form: 


Sij  "  '  P6ij  +  Tij 


(11) 


The  dependent  variables  in  the  conservation  equations  can  be  ex¬ 
panded  into  the  following  form: 


u  »  u  +  u'  (12a) 
v  ■  v  +  v'  (12b) 
P  »  P  +  P'  (12c) 
P  ■  P  +  p'  (12d) 
e  ■  e  +  e*  ( 1 2e ) 
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In  these  expansions  the  barred  variables  represent  time  averaging  over 
a  time  interval  that  is  long  compared  to  turbulent  eddy  fluctuations,  yet 
small  compared  to  macroscopic  flow  changes.  The  primed  variables  then 
represent  fluctuations  due  to  the  turbulent  nature  of  the  flow.  As  dis¬ 
cussed  by  Chapman  (Reference  18),  this  time  averaging  approach  is  valid 
since  the  frequencies  of  most  unsteady  flows  of  interest  are  a  factor  of 
10  to  100  below  the  mean  frequency  of  the  turbulent  eddy  motion. 

If  the  dependent  variables  u,  v,  and  e  are  mass  averaged  as 
described  in  Reference  19,  and  p  and  P  are  mean  (time  averaged)  state 
variables,  then  the  conservation  equations  can  be  written  in  the  form  of 
mean  flow  equations  as: 

P.t  +  (PuJ.j  -  0 

<pSi>.t  +  KPVj)  +  It  -  (Tjj  -  pujuj)],^  -  0 
(pe),t  +  [peilj  +  +  puje'  -  (t^  -  pu’uj) ] ,  -  0 

where  a  higher  order  mean  energy  dissipation  term  in  u.!  has  been  neglected 
in  the  energy  Equation  15. 

The  term  [-pujuj]  is  known  as  the  Reynolds  stress.  It  represents 
a  momentum  transfer  caused  by  turbulent  fluctuations  present  in  the  flow- 
field.  This  Reynolds  stress  term  can  be  written  as  an  apparent  stress 
caused  by  the  turbulent  nature  of  the  flow: 


(13) 

(14) 

(15) 


ij 'turb 


Puiuj 


(16) 


Since  air  is  essentially  an  isotropic  fluid,  the  mean  stress  term  can 
be  expanded  into  its  normal  and  shear  stress  components  as: 


T«  '  Juk.k  ‘ij + "  <“i.j  *  “j ,,> 


(17) 
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The  turbulent  stress  term  can  then  be  written  in  analogous  form  as: 


Tij  ^turb 


Uk,k  6ij 


(18) 


where  and  e  are  the  turbulent  viscosity  coefficients  of  the  flow. 

The  coefficient  c  is  known  as  the  eddy  viscosity,  and  is  analogous  to 
the  molecular  viscosity  coefficient  p.  However,  e  is  more  a  property  of 
the  dynamics  of  the  flow,  whereas  p  is  only  a  property  of  the  fluid. 
Combining  the  mean  and  turbulent  stress  terms,  an  overall  stress  term 
can  be  written  as: 


total 


(X+V  \,k  6ij +  (y+£)  (Gi,j  +  =jfi> 


(19) 


In  the  energy  equation  an  additional  unsteady  term  appears.  This 
term  is  by  nature  an  apparent  heat  flux  caused  by  the  fluctuations 
inherent  to  turbulent  flow  and  can  be  written  as: 


,  *=  pu'e' 
turb  j 


(20) 


If  the  heat  flux  term  is  defined  by  the  Fourier  heat  equation  as: 


qj  =  "kT,j  =  "  (cpu/Pr) 


(21) 


then  by  the  former  analogy  qj  tur^  can  be  written  as: 


turb  '  -  «WV  ’'j 


(22) 


where  again  l  is  the  eddy  viscosity  coefficient,  and  Pr^  is  the  turbulent 
Prandtl  number  of  the  flow.  Combining  these  two  heat  fluxes,  a  total 
heat  flux  can  be  written  as: 


qJ  total 


-Cp  (w/Pr  ♦  t/Prt)  Tfj 


(23) 
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The  mean  conservation  equations  can  then  be  written  in  the  follow¬ 
ing  form,  where  the  overbars  on  the  terms  are  dropped  for  convenience,  and 
where  the  values  cf  the  shear  stresses  and  heat  fluxes  are  the  total 
values : 

P.t  +  “  0  (24) 

(pUi).t  +  [(pui)uj  +  -  o  (25) 

(pe),t  +  [(pe)  uj  +  ^  Ti^*j  “  0  (26) 

Since  the  flowfields  of  interest  are  either  two-dimensional  or 
axisymmetric  in  nature,  the  mean  conservation  equations  can  be  written 
in  the  following  compact  vector  form: 

i£  +  il  +  jL_lkl!Gi  _  .  jl  (27) 

3t  3X  +  rj„  3r  rje  v  1 

where  jQ  =  0  or  1  for  either  the  two  dimensional  or  axisymmetric  cases, 
respectively,  and 


PU2 

pu  -  0 

»  F  ■ 

XX 

PUV  -  T 

pue  +  $  -  uo  -  vt 

L  x  xx  xr 

_  PUV  -  T 

G  *  2  xr 

pvz  -  a 

pve  +  q  -  ut  -  vo 

r  xr  rr, 


where 


-P  +  (X+X  )  div  u  +  2(u  +  e)  — 


°rr  -  -P  +  (X+Xt)  div  u  +  2 (u  +  t)  Y- 
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oH  -  -P  +  (A+Xt)  div  u  +  2  (p  +  e)  ^ 
txr-  (u  +  c)  (£  +  £) 


•  _  r  fV_  ,  i _ s 

qx  P  '■Pr  Prt;  3x 

•  _  .p  £ .  3T 

qr  “  -Cp  (Pr  Pr^  3r 


-*•  3u  ,  3v  .  .  v 

div  u«  —  +  — +  j„  — 
3x  3r  J  r 


The  coefficient  of  viscosity  u  for  air  can  be  assumed  to  vary 
according  to  Sutherland's  law  (Reference  19): 

p  -  (2.27  •  10-8)T1,5/(T  +  198.6)  (lbf  -  sec/ft4) 


The  second  coefficient  of  viscosity  is  assumed  as  the  following 
by  applying  Stokes  Hypothesis: 

(A  +  A£)  -  -2/3  (p  +  e)  (3?) 

The  governing  equations  for  the  problem  of  interest  now  consist 
of  the  four  conservation  equations  in  matrix  form  (Equation  27)  with 
four  unknown  dependent  variables  p,  pu,  pv  and  pe.  The  perfect  gas 
law  is  used  to  define  the  pressure  in  terms  of  these  conservative 
variables,  and  a  model  of  the  dependence  of  the  eddy  viscosity  on  the 
mean  flow  must  be  introduced  to  overcome  the  "turbulent  closure"  problem 
inherent  in  the  turbulent  mean  conservation  equations. 

For  numerical  computation  in  a  transformed  (f.,n)  Cartesian  plane, 
the  matrix  form  of  the  conservation  equations  (Equation  27)  can  be 
written  as: 


♦  "s.  Is  ♦ -b  \  ^  * 
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where  £  and  n  are  now  the  independent  variables,  and  the  transformation 
derivatives  £  ,  Cr,  nx»  and  or  are  obtained  numerically  from  a  mapping 
procedure.  Equation  38  is  actually  in  weak  conservation  form  due  to 
the  varying  coefficients  in  front  of  the  derivatives,  and  also  due  to 
the  source  term  in  the  axi symmetric  case. 

2.  BOUNDARY  AND  INITIAL  CONDITIONS 

Boundary  and  initial  conditions  must  be  given  in  order  to  solve  the 
conservation  equations  which  govern  the  flowfield.  These  conditions  must 
be  carefully  specified,  since  many  flow  features  such  as  shock  waves, 
boundary  layers,  and  recirculation  areas  arise  from  boundary  conditions. 
For  the  solution  of  a  symmetric  two-dimensional  or  axi symmetric  supersonic 
jet  embedded  in  a  supersonic  external  flowfield,  the  domain  of  interest 
can  be  defined  as  shown  in  Figure  6.  Only  one-half  of  the  total  nozzle 
flowfield  needs  to  be  considered  due  to  the  axis  of  symmetry  on  the  jet 
centerline.  The  remainder  of  this  chapter  will  detail  the  specific 
boundary  conditions  that  are  pertinent  to  this  problem. 


Figure  6.  Physical  Domain  for  the  Computational  Solutions 
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a.  The  Upstream  Boundary 

Inflow  conditions  on  the  upstream  jet  boundary  (AB)  and  the 
upstream  external  boundary  in  the  freestream  (CD)  are  completely  specified. 
Velocity,  pressure,  and  temperature  profiles  determined  from  auxiliary 
computations  or  known  experimentally  on  this  boundary  fix  p,  pu,  pv, 
and  pe  for  the  duration  of  the  problem  solution  at  the  boundary. 

b.  The  Upper  Boundary 

The  upper  boundary  (DE)  must  accurately  represent  a  free  flight 
condition  where  mass  flow  is  allowed  across  this  boundary  embedded  in 
the  supersonic  external  flowfield.  Weak  shock  waves  and  Prandtl -Meyer 
expansion  waves  must  also  exit  this  boundary  without  being  reflected 
artificially  back  into  the  domain  of  interest.  One  condition  that 
allows  this  is  the  assumption  of  a  simple  wave  solution: 


-  0 

DE 


(39) 


where  Xj  is  the  straight  left  running  characteri stic  line  passing  through 
each  point  on  the  upper  boundary.  This  characteristic  line  is  determined 
only  by  the  value  of  the  Mach  number  and  flow  angle  of  the  supersonic 
flow  present  near  the  upper  boundary.  This  condition  assumes  that  the 
flow  along  this  boundary  is  inviscid  and  homentropic. 


c.  The  Downstream  Boundary 

The  downstream  boundary  (EF)  is  unique  in  that  no  rigorous 
assumptions  can  be  made  about  either  the  variables  or  their  gradients 
unless  the  boundary  is  placed  a  great  distance  downstream.  In  this  case 
a  no  change  condition 


9f 

dx 


(40) 


could  be  assumed,  where  f  denotes  the  primitive  variables  p,  u,  v,  and 
T.  For  the  case  where  the  downstream  boundary  is  placed  where  gradients 
do  exist,  an  extrapolation  method  based  on  this  fact  can  be  reasonably 
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applied.  One  such  method  is  to  assume  that  a  flow  gradient  accurate  to 
second  order  can  exist.  This  can  be  implemented  as: 


.  333f. 

(Ax  - r-) 

3x 


EF 


(41) 


In  other  words,  gradients  can  occur  which  are  parabolic  with  respect  to 
x.  This  condition  is  reasonable  if  the  gradients  at  this  boundary  are 
not  severe  as  in  the  case  where  a  strong  shock  wave  exits  the  boundary. 

d.  The  Centerline 

The  centerline  boundary  (AF)  is  a  line  of  symmetry  with  no  mass 
or  energy  flux  across  it.  Therefore,  the  following  boundary  conditions 
can  be  appl ied 


V 


0 


(42) 


3p 

3r 


(43) 


3u 

3r 


(44) 


Since  the  v  component  of  velocity  is  zero  on  the  centerline,  this 
boundary  is  also  a  streamline  in  the  jet  flow.  For  steady,  adiabatic 
flow  with  negligible  volume  forces,  the  total  enthalpy  along  any 
streamline  is  a  constant.  Therefore,  along  the  centerline, 

■  constant  (45) 


AF  "  (CPT  +  ¥ 


Since  the  condition  at  the  jet  exit  is  specified,  the  centerline  boundary 
can  be  properly  posed  using  this  approach. 
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e.  Nozzle  Walls 

The  nozzle  walls  (BG,  GH,  and  CH)  are 
impermeable  boundaries.  This  assumption  gives 


considered  to  be  no-slip, 
the  conditions  that: 


u|wall  ®  ( 46a ) 

and 

v Iwall  “  0  (4ob) 


Since  the  stainless  steel  nozzle  consists  of  thin-walled  material 
with  a  thermal  conductivity  much  greater  than  that  of  the  surrounding 
fluid,  the  nozzle  walls  are  assumed  to  be  at  a  constant  temperature: 


.  (47) 

T  -  constant 

|wall 

This  wall  temperature  is  determined  by  applying  a  heat  flux 
balance  across  the  jet  and  freestream  boundary  layers  as  outlined  in 
Appendix  A. 

The  pressure  on  each  nozzle  wall  is  unknown,  but  can  be 
approximated  by  applying  the  degenerate  form  of  the  appropriate  normal 
momentum  equation  at  each  nozzle  surface  to  obtain  the  following: 


3P 

3n 


wall 


3t 


sn 

3s 


wall 


0 


(4b) 


In  this  expression  n  is  the  directional  normal  to  the  wall  surface,  and 
s  is  the  direction  parallel  to  the  surface. 


f.  Initial  Conditions 

Since  the  governing  equations  contain  time  dependent  terms, 
initial  conditions  must  be  specified  before  the  solution  process  can  begin. 
The  specification  of  these  initial  conditions  is  somewhat  arbitrary  in 
nature,  although  steep  gradients  must  be  avoided  to  prevent  numerical 
divergence  during  the  solution  process.  Since  the  flow  is  predominantly 
supersonic  in  nature,  the  incoming  flow  profiles  will  have  a  dominant 
effect  on  the  solution  in  the  whole  computational  domain.  The  incoming 
flow  conditions  are  imposed  on  the  complete  domain  as  discussed  in  the 
section  on  initial  condition  implementation  of  Section  IV. 
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SECTION  III 
NUMERICAL  PROCEDURE 

The  numerical  procedure  consists  of  solving  the  governing  equations 
with  applicable  boundary  and  initial  conditions  through  the  use  of 
appropriate  finite  difference  techniques  on  a  high-speed  computer.  This 
procedure  can  be  broken  down  into  several  elements  which  include  the 
finite  difference  coordinate  system,  the  solution  algorithm,  and  the 
convergence  criteria  used  in  the  computational  solution.  Each  of  these 
areas  will  be  discussed  in  this  chapter. 

1.  COORDINATE  SYSTEM 

a.  Domain  of  Computation 

The  physical  domain  of  computation  consists  of  a  rectangular  area 
defined  by  orthogonal  coordinates  (x,r)  as  shown  in  Figure  7.  The  mesh 
consists  of  IL  points  on  the  x  axis  and  JL  points  on  the  r  axis,  where 
IL  and  JL  are  dependent  on  the  extent  of  the  physical  domain  required 
for  the  particular  jet  plume  case  of  interest. 

b.  Mesh  Stretching 

In  order  to  obtain  an  accurate  numerical  solution  of  a  viscous 
flowfield,  the  mesh  spacing  must  be  much  finer  in  areas  containing 
relatively  high  gradients  of  the  variable  properties  such  as  velocity, 
density,  and  temperature.  In  the  nozzle  flowfield  these  high  gradient 
areas  include  the  boundary  layers  on  the  nozzle  walls  and  the  shear 
layer  in  the  wake  of  the  nozzle  annulus.  This  stretching  is  accomplished 
through  the  use  of  a  patched  exponential  stretching  scheme  of  the  follow¬ 
ing  form: 

Cn(j) 

r(j)  *  L  (^-r - -1)  for  J  -  1,N  (49) 

<e  -  1) 


where  L  and  n  are  defined  as  shown  in  Figure  8.  The  constant  C  is  deter¬ 
mined  by  the  minimum  spacing  Ar  .  desired  for  the  mesh  next  to  the  wall 
boundary.  Applying  the  desired  Ar  .j  to  Equation  49  gives 


C 


1_ 

An 


Loge 


(1  + 


Ar 


min 


L 

m 


(eC  -  1» 


(50) 
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The  value  of  C  is  then  obtained  through  the  use  of  an  iterative  Aitken 
extrapolation  technique. 

This  mesh  stretching  procedure  is  applied  in  the  radial  direction 
on  both  sides  of  the  nozzle  wall  where  boundary  layers  are  present.  It  is 
also  applied  in  the  axial  direction  at  the  end  of  the  nozzle  where  the 
jet  flow  begins  to  expand  or  contract  and  where  the  near  wake  due  to  the 
nozzle  annulus  begins  to  form. 

c.  Adaptive  Mesh 

It  is  desireable  that  the  fine  mesh  remain  in  the  areas  of 
relatively  high  velocity  and  temperature  gradients  as  the  solution  pro¬ 
gresses  towards  convergence.  This  is  not  a  problem  in  the  case  of 
boundary  layers  that  are  adjacent  to  a  fixed  wall,  but  is  a  concern  in 
the  free  shear  layer  area  generated  by  the  nozzle  wake  and  the  interaction 
between  the  jet  and  freestream  flows.  This  shear  region  on  the  jet  plume 
boundary  will  deflect  to  a  degree  that  is  primarily  dependent  on  the 
nozzle  pressure  ratio.  The  fine  mesh  region  should  therefore  also  be 
adapted  to  conform  to  the  general  position  of  the  shear  region. 

Hirt  (Reference  21)  has  used  a  technique  in  the  solution  of  free 
surface  flows  that  allows  the  grid  to  adapt  as  the  solution  progresses. 

The  following  kinematic  equation  is  applied  in  the  region  where  the  shear 
layer  is  present. 


3r 

3t 


CA  (V  -  U 


(51) 


This  equation  insures  the  condition  that  as  the  solution  con¬ 
verges,  the  physical  slope  of  the  constant  n  finite  difference  cell 
boundaries  is  the  same  as  that  of  the  velocity  vectors  near  each  cell. 

When  applied  in  a  finite  difference  format,  the  grid  can  then  adapt  to 
the  placement  of  the  shear  layer  as  shown  in  Figure  9.  Details  of  this 
process  are  explained  in  Appendix  B. 

d.  Coordinate  Transformation 

The  physical  domain  as  typically  shown  in  Figure  9  is  mapped  to  a 
unit  square  in  the  computational  plane  shown  in  Figure  10.  The  constant 
r.  lines  are  aligned  parallel  to  the  centerline  and  the  constant  t  lines 
are  in  the  direction  normal  to  the  centerline.  The  numerical  algorithm 
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operates  on  this  coordinate  system  using  the  transformed  conservation 
equations  (Equation  38).  Care  must  be  taken  in  generating  the  physical 
mesh  so  that  smoothness  of  the  transformation  coefficients  (cx,  £r,  n 
and  'i  )  is  retained  in  order  to  reduce  numerical  errors  caused  by  the 
mesh  configuration. 

2.  SOLUTION  ALGORITHM 

a.  MacCormack's  Method 

The  weak  conservative  form  of  the  two-dimensional,  time-dependent 
Navier-Stokes  Equations  (Equation  38)  is  solved  using  MacCormack's  ex¬ 
plicit  finite  difference  method  (Reference  22).  This  algorithm  is  an 
efficient  Lax-Wendroff  type  differencing  scheme  of  second  order  accuracy 
which  utilizes  time-splitting  and  two  step  predictor-corrector  techniques. 
MacCormack's  algorithm  was  chosen  for  application  to  the  nozzle  problem 
because  of  its  previous  success  in  computing  inviscid-viscous  interacting 
flows,  its  stability  in  supersonic  flow,  and  its  computational  efficiency 
achieved  by  time-splitting  the  finite  difference  operators. 

The  computational  solution  is  advanced  in  time  by  applying  the 
numerical  operator  to  the  solution  of  the  flowfield  at  time  t.  This  can 
be  written  as: 


ua,  n,  t  +  At)  -  L(At)  •  U(£,  n,  t) 

where  L(&t)  is  the  two-dimensional  numerical  operator  representing 
MacCormack’s  algorithm  acting  on  the  transformed  conservation  equations. 
Through  the  use  of  a  time-splitting  technique,  this  two-dimensional 
operator  L(At)  is  separated  into  two  one-dimensional  sweep  operators  in 
the  and  n  directions.  The  operator  Lr(At)  denotes  the  solution  of 
the  equation: 


3U  +  ,  JF  S(rjeG 

3t  *  Sar^,  S  H 


0 


(53) 
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in  the  £  direction  by  time  increment  At  seconds.  Similarly,  the 
operator  L^At)  represents  the  solution  of 


3U 

at 


+ 


a(rj6G) 

3n 


(54) 


in  the  n  direction  by  a  time  increment  At  seconds.  The  dependent 
variable  vector  U(c,  n»  t)  can  then  be  advanced  in  time  as 


U(£,  n,  t  +  At)  -  [L^M/2(At/M)-Ln(At)-L5M/2(At/M)]-U(C,  n,  t)  (55) 
wi  th 


At  -  At.  if  At.  <  At 
C  C  n 


or  as 


U(e.  n.  t  ♦  At)  -  [LnN/2(At/N)*L£(At)-LnN/2(At/N)].U  (£,  n,  t)  (56) 


with 


At  -  At  if  At  <At_ 
n  H  C 


In  these  equatioi  M  and  N  are  the  smallest  even  integers  of  the  quotients 
(At^/At^)  and  (a  'At^),  respectively,  and  At.  and  At^  are  the  maximum 
allowable  time  s'  ps  in  the  £  and  n  directions  as  determined  by  the 
Courant-Friedrichs-Lewy  (CFL)  limit  discussed  in  the  next  section  on 
stability.  The  values  of  M  and  N  are  usually  equal  to  two  for  the  grid 
distribution  used  in  the  solution  of  the  nozzle  problem.  This  sets  up 
a  truly  alternating  direction  procedure  that  is  desirable  when  gradients 
exist  in  more  than  one  direction. 

The  finite  difference  forms  of  the  sweep  operators  consist  of  a 
predictor-corrector  procedure  which  increases  the  accuracy  of  the  time- 
dependent  term  evaluations.  This  method  utilizes  one-sided  differencing 
in  the  direction  of  sweep,  but  central  differencing  in  the  direction 
perpendicular  to  the  sweeping  coordinate.  At  the  completion  of  the 
predictor-corrector  process,  this  method  is  equivalent  to  a  second  order 
central  differencing  scheme  in  two  dimensions. 
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The  Lr  (At)  sweep  operator  represents  the  following  numerical 
procedure.  A  predicted  intermediate  value  is  computed  by  the  expression: 


u"  -  77  1(0,  ,  (f"  , 

i,j  x  l,j 


F 


n 

i-l.j 


) 


rl' 
i.j 


‘Vid  ‘'id  Gli 


-  r 


i-l.j  i-l.j 


)] 


(57) 


where  ll1?  .  is  the  known  value  at  time  t,  and  is  the  intermediate 

1  *  J  '  y  J 

predictor  value.  The  actual  computed  value  at  time  (t+At)  is  then 
calculated  by  applying  the  following  corrector  algorithm: 


"tj  '  1/2  <  "id  +  -  <f>  1‘Vi.j  -  >£?> 


r j o  rn+H 

r--  ’r'i.j  '  i+l,j  i+l.j 
i.j 


1 

_J  o 


+  <U,  .  <r H,  ,  G‘17  ,  -  gTDU 


i.j  i.J 


(58) 


In  this  c  sweep  predictor-corrector  algorithm,  the  matrices  F 
and  G  are  functions  of  the  following  difference  quotients: 


-  U?  .  U? 


-  U, 


F®  Qm  m  £  (u®  ...i+.Lt J — _ Ll  i > _ i .J-l) 

i.j  ’  i.j  6  ^  i.j’  AC  ’  2An  } 


(59) 


This  scheme  can  be  depicted  graphically  as  shown  in  Figure  11. 

The  (At)  numerical  sweep  operator  is  formulated  in  a  similar 
manner.  The  intermediate  predictor  value  is  given  by  the  expression: 


<J  -“VO  "Vi.)  "id  -  Fid-1> 


<Vi.i  "id  ci.i  -  rid-i  'id-i” + 


jo  At  H 


LI 


(60) 


i.j 
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£  SWEEP 

PREDICTOR  CORRECTOR 

j 

j-1 


r\  SWEEP 

PREDICTOR  CORRECTOR 

j  +  1 

j 

j-1 

i-l  i  i +1  i -1  '  * 

Figure  11.  Graphical  Representation  of  the  Numerical  Sweep 
Operators 
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The  corrector  value  at  Lime  (t+At)  is  then  given  by: 


u”+;  -  1/2  {  u”  +  uf '1  -  <|£->  Kn  ),  .  -  f"+?) 

i.j  i,j  i.j  An  x  i ,j  i.J+l  i,3 

■1  A  r  Hn+^ 

+  JL.  (n  \  J»  cn+!*  _  Jo  Gn+J*n  +  3°  1 J 

Jo  Vi,j  (ri,j+l  Gi,j+1  ri,j  Gi,j)]  +  1. 

i.j 


The  matrices  F  and  G  are  now  functions  of  the  following  difference 
quotients : 


(61) 


_  U  _  ti”1 

•r  cm '  K.r  ' 1-1  ,J  •  ‘'J> 


(62) 


This  scheme  is  also  depicted  graphically  in  Figure  11. 

The  MacCormack  algorithm  satisfies  the  requirements  of  stability, 
consistency,  and  accuracy  if  the  following  conditions  are  met:  It  is 
stable  if  the  time  step  limit  satisfies  the  CFL  stability  condition 
given  in  the  next  section.  It  is  consistent  if  the  sum  of  the  time 
steps  for  each  the  sweep  operators  is  equal,  and  the  algorithm  is 

second  order  accurate  if  the  sequence  of  sweep  operators  is  symmetric. 

b.  Stability  Condition 

Although  the  stability  of  the  finite  difference  equations  cannot 
be  analyzed  when  applied  to  the  complete  set  of  Navier-Stokes  equations, 
MacCormack  and  Baldwin  (Reference  23)  have  shown  that  considerable  insight 
can  be  gained  by  separately  analyzing  the  inviscid,  diffusion,  and  mixed 
derivative  parts  of  these  equations.  Stability  criteria  are  determined 
by  calculating  the  eigenvalues  of  the  associated  ampl i f ication  matrix 
for  each  part.  If  this  eigenvalue  procedure  is  applied  to  the  inviscid 
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terms  of  the  Navier-Stokes  equations  for  flow  in  a  Cartesian  coordinate 
system,  the  following  stability  conditions  emerge: 


(63) 

(64) 


where  c  is  the  local  speed  of  sound.  Consideration  of  only  the  diffusive 

2  2  2  2 

terms  which  contain  a  U/ax  and  a  U/ar  in  the  complete  Navier-Stokes 
equations  gives  the  conditions: 


Atx  1  1/2 


X 
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(ii  +  X  ) 

rr  Pr  ' 


(65) 


and 


Atr  <  1/2 


*<£  + 

p  Pr 


(66) 


Finally,  analysis  of  the  mixed  derivative  terms  found  in  the  Navier- 
Stokes  equations  gives  the  condition: 


At  -  At  <  Ax  Ar 

x  r  -  - - 

j  [-(X  +  A  >  <u  +  t)P 


(67) 


For  the  finite  difference  equations  applied  to  the  complete  set  of 
Navier-Stokes  equations,  the  stability  criteria  can  then  be  estimated  as: 

At  <_  minimum  (At  ,  At  )  (58) 

i.J 
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where 

At  < 
x  — 


Ax 


M  +  '  ♦  h {  S  <p?  +  fet>  +  h  ><"  +  Ol^j 


(69) 


and 


At  < 


Ar 


r  m + ' + 1 { <jf  +  fet> +  k  +  or) 


(70) 


For  the  present  non-Cartesian  jet  flow  cases  that  were  computed, 
the  maximum  time  step  was  calculated  as: 


At  -  K  .minimum  (At,,  At  ) 
cfl  5  " 


(71) 


where 
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(73) 


where 


as4-  ^i.j 

*  "i-l,/  +  <‘l.J  - 

(74) 

ASn" 

-  W*  *  "l.3  -  W*1' 

(75) 
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(ui.j'  vi,j'  "i-l,r 
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u  “  Maximum 

<ui,j-  V1,J’  "i.J-l- 

(77) 
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CFL  factors,  K  ,  varying  between  0.35  and  0.80  were  used  during  the 
computations.  Flow  solutions  involving  relatively  large  viscous  wakes 
containing  recirculation  regions  required  much  smaller  allowable  time- 
steps  with  CFL  factors  on  the  order  of  0.35-0.40. 


c.  Numerical  Damping 

Strong  shocks  imbedded  in  a  flowfield  being  solved  computationally 
can  often  cause  numerical  oscillations  (Reference  17)  which  may  lead  to 
program  failure  due  to  physically  unrealistic  values  of  computed  pressure, 
density  or  temperature  (Gibb's  phenomenon).  These  oscillations  are 
caused  by  numerical  truncation  errors  and  can  be  reduced  by  refining  the 
grid  in  the  areas  of  shock  locations.  However,  this  can  be  impractical 
when  the  oscillations  are  of  a  transient  nature  caused  by  computational 
start-up  or  re-start  procedures,  or  where  the  shock  location  varies  for 
different  experimental  cases  and  mesh  refinement  for  each  individual  case 
is  undesirable.  In  this  situation,  a  fourth  order  pressure-gradient 
damping  concept  as  introduced  by  MacCormack  and  Baldwin  (Reference  23) 
can  be  applied  to  increase  the  stability  of  the  numerical  algorithm  by 
adding  artificial  viscosity. 


This  damping  scheme  is  applied  in  both  the  i  and  n  directional 
sweeps.  In  the  numerical  sweep,  damping  is  implemented  by  replacing 

+  GQ  )  in  Equation  53. 


Fii,j  (Fii,0  +  F0U  >  and  Gii,j  ^  '  Gf  i .  j 

The  predictor  and  corrector  steps  in  this  case  are  represented  by 
ii=i  and  ii  =  i+1,  respectively.  The  damping  terms  F^  and  Gp  are  in 
the  following  form: 


ii.J 
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and 


(Ui+lt 3  "  Ui.J) 


'ii.J 


a  ( |v  I  +  c  )  l-LifJ+l"  2Pli,J  *  Pii,J-l. 
1,1  Ai'j  U’J  <Pii.J+l  +  2Pii.3  + 
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(79) 
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In  the  n  sweep,  a  similar  procedure  is  implemented  by  replacing  F.  .. 

*  »JJ 

)  where 


by  (Fi,Jj  +  FDi  .  >  a"d  Gi,jj  by  tGi,jj  +  GD,  „ 
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(81) 


In  this  case  the  predictor  and  corrector  steps  are  represented  by  jj  =  j 
and  jj  =  j+1 ,  respectively.  For  both  sweeps  a ^  and  are  damping 
constants  where  norma lly 


0.5 


(82) 


This  damping  technique  produces  fourth  order  viscosity  like 
terms  of  the  form 


|ui+c 
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which  are  added  to  the  difference  equations.  The  magnitude  of  these 
damping  terms  is  proportional  to  the  second  derivative  of  pressure  and 
is  significant  only  in  regions  of  pressure  oscillation  where  the  trun¬ 
cation  error  is  already  adversely  affecting  the  solution. 

3.  CONVERGENCE  CRITERIA 

Convergence,  as  applied  in  this  section,  refers  to  iteration  con¬ 
vergence  as  opposed  to  truncation  convergence,  which  involves  the 
convergence  of  the  solution  of  the  FOE  to  the  solution  of  the  PDE  as 
ax,  Ar,  and  At  -*•  0.  Iteration  convergence  refers  to  the  arrival  at  a 
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solution  to  the  discretized  Navier-Stokes  equations  within  some  acceptable 
tolerance  through  the  use  of  an  iterative  process.  As  stated  by  Roache 
(Reference  17),  there  are  no  definitive  criteria  for  iteration  convergence. 

A  somewhat  subjective  judgement  of  convergence  must  be  made  based  upon  an 
examination  of  the  iterative  behavior  of  the  solution  flow  variables. 
Different  flow  variables,  as  well  as  variables  at  different  locations, 
converge  at  different  rates.  If  the  slowest  converging  variable  in  the 
flowfield  is  known,  it  should  be  the  most  closely  examined  for  convergence. 

In  the  present  case  for  a  coflowing  supersonic  nozzle  with  a 
relatively  thick  base  annulus,  an  examination  of  the  flow  variables  revealed 
that  the  slowest  converging  variable  was  the  base  pressure  of  the  nozzle 
annulus.  The  location  of  this  base  pressure  is  within  the  subsonic  flow 
area  involving  recirculation  in  the  wake  of  the  nozzle  annulus  as  shown 
in  Figure  12.  The  flow  variables  in  this  subsonic  region  converged  much 
more  slowly  than  did  those  in  the  predominantly  supersonic  jet  and  free- 
stream  flows. 

Since  the  coflowing  nozzle  problem  primarily  involves  high  Reynolds 
number  flow,  the  advective  terms  in  the  conservation  equations  dominate 
the  viscous  diffusion  terms.  A  characteristic  time  for  a  disturbance  to 
cross  the  flowfield  may  then  be  characterized  by: 

Cch  u  (83) 

ch 

where  L  is  the  length  of  the  flowfield  in  the  direction  parallel  to  the 
characteristic  velocity  uc^.  For  the  jet  problem  uc^  was  represented 
by  um.  Since  in  general  the  magnitude  of  uu  was  less  than  u^,  this  gave 
a  more  conservative  estimate  of  the  characteristic  time. 

The  convergence  criteria  was  then  established  by  the  following 
procedure.  The  numerical  solution  was  either  initially  started,  or  re¬ 
started  from  a  previous  case.  As  the  solution  converged,  the  base  pres¬ 
sure  was  monitored  until  its  magnitude  varied  less  than  +  IT  for  one 
characteristic  time  period.  At  the  end  of  this  characteristic  time  period 
the  solution  was  stopped  as  convergence  was  achieved.  Visual  comparison 
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of  Mach  Number  and  density  profiles  over  the  flowfield  confirmed  the 
convergence  of  the  solution  using  this  procedure.  A  sample  base 
pressure  convergence  plot  is  shown  in  Figure  13. 


Figure  12.  Base  Pressure  Tap  Location  in  the  Flowfield 
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SE(T,ON  IV 

BOUNDARY  AND  INITIAL  CONDITION  IMPLEMENTATION 

As  in  Section  II,  boundary  and  initial  conditions  must  be  defined 
in  order  to  solve  the  conservation  equations  which  govern  the  flowfield. 
Values  of  the  dependent  variables  for  points  on  the  boundaries  of  the 
computational  domain  must  be  specified  in  order  for  the  numerical 
operators  to  solve  the  flowfield  correctly.  This  section  presents  the 
numerical  formulation  of  the  boundary  and  initial  conditions  used  for 
the  solution  of  the  axisymmetric  nozzle.  The  conditions  were  presented 
in  a  mathematical  context  in  Section  II. 

1.  THE  UPSTREAM  BOUNDARY 

The  flow  properties  on  the  upstream  boundary  (AB  and  CD  of  Figure  6) 
are  held  fixed  for  the  duration  of  the  computational  solution.  The  values 
of  these  properties  were  derived  in  the  following  manner.  In  the  external 
flowfield,  a  parabolized  Navier-Stokes  solution  (Reference  24)  was  com¬ 
puted  for  the  ogive  body  used  in  the  experimental  nozzle  tests  as  shown 
in  Figure  14.  This  solution  determined  that  the  pressure  gradient  at  the 
inflow  boundary  in  the  external  flow  stream  is  negligible,  and  that  the 
static  pressure  at  the  inflow  boundary  is  99%  of  that  in  the  undisturbed 
flow  in  the  wind  tunnel.  The  static  pressure  along  the  ogive  body 
surface  shown  in  Figure  15  was  then  used  as  an  input  to  a  two-dimensional 
turbulent  boundary  layer  code  (References  25  and  26)  along  with  the  other 
specified  freestream  conditions  (M^  T0co,  Re^,  T  )  to  produce  input  velocity 
component  and  temperature  profiles  at  the  upstream  boundary.  The  con¬ 
servative  variables  on  this  boundary  (CD)  are  then  calculated  using  these 
profiles  ar.d  the  static  pressure  along  the  boundary.  Flow  variables  on 
the  upstream  boundary  in  the  jet  flow  (AB)  are  determined  in  a  similar 
manner.  The  same  boundary  layer  code  is  applied  using  the  jet  exit 
conditions  and  the  length  from  the  nozzle  throat  to  the  nozzle  exit 
plane  as  the  boundary  layer  starting  length.  Again,  profiles  for  the 
velocity  components  and  temperature  are  obtained  along  the  boundary. 

Values  for  the  conservative  variables  are  then  computed  using  these 
profiles  and  the  value  of  the  pressure  at  the  jet  exit.  Since  the  value 
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of  the  vertical  velocity  component  is  zero  on  the  centerline  boundary,  a 
polynomial  fit  is  used  to  set  the  vertical  velocity  profile  from  the  edge 
of  the  boundary  layer  to  the  centerline. 

2.  THE  UPPER  BOUNDARY 

The  upper  boundary,  labeled  DE  on  Figure  6,  utilizes  the  simple  wave 
procedure  outlined  by  Roache  {Reference  17,  pp  282-283).  This  procedure 
assumes  that  properties  are  constant  along  a  straight,  left-running 
characteristic  line  passing  through  each  point  on  the  upper  boundary. 

The  position  of  this  line  running  through  a  boundary  point  (i,JL)  is 
determined  by  the  angle  UM  +  0),  where 

^  ■  arcsin  (1/M^  (84) 

is  the  local  Mach  angle  for  supersonic  flow,  and 

0  -  arctan  (v/u)  (gg) 

is  the  local  flow  direction.  The  properties  on  this  characteristic 
line  are  determined  by  linear  interpolation  involving  the  properties  at 
points  (1-1 ,  JL),  (i-1,  JL-1),  and  (i ,  JL-1 )  as  shown  in  Figure  16. 

Points  (i-1,  JL-1)  and  (i,  JL-1)  are  points  interior  to  the  computational 
domain,  and  point  (i-1,  JL)  is  either  a  point  on  the  inflow  boundary 
(i  =  1)  or  a  previously  defined  point  on  the  upper  boundary  resulting 
from  the  left  to  right  sweep  along  the  boundary  using  this  technique. 

As  shown  in  Figure  16,  one  of  two  possible  interpolation  schemes  is 
applied  depending  on  the  local  values  of  the  quantities  (uM  +  o)  and 
Ur/ ax).  For  the  case  of  tan  UM  +  o)i_1  JL_1  >  Ar/Ax,  the  position 
of  the  characteristic  line  lies  between  the  points  (i-1,  JL-1)  and 
(i,  JL-1).  The  properties  at  the  point  p,  and  thus  at  point  (i,  JL), 
can  then  be  determined  by: 

l 

fi,JL  fp  “  f i-1 , JL-1  +  ^ (fi,JL-l  '  fi-l,JL-l)  (86) 
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Figure  16.  Upper  Boundary  Condition  Schematic 
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The  value  of  a,  and  thus  the  position  of  the  point  p,  can  be  determined 
by  the  following  procedure.  If  the  quantity  w  is  defined  by: 

w  ■  tan  [90°  -  (l^  +  ©)]  (87) 


then  by  geometry: 

w  -  (Ax  -  £)/ Ar 
P 


(88) 


If  the  interpolation  procedure  of  Equation  86  is  applied,  then 


wp  “  +  (AxJ  (wi,JL-l  ~  Wi-1,JL-1) 


(89) 


Equating  Equation  88  to  _quation  89  and  solving  for  i  gives: 
(Ax/Ar)  -  v. 


I  m 


(w 


i-1, JL-1 


i.JL-l  "  Wi-1,JL-1)/Ax  + 


(90) 


For  the  alternate  case  of  tan  (um  +  o)^_^  JL_-|  <  (Ar/Ax),  the 
position  of  the  characteristic  line  lies  between  the  points  (i-1,  JL-1 ) 
and  (i-1,  JL).  The  properties  at  p  are  then  determined  by  the  equation: 


i,  JL 


“  fp  “  f i-1, JL-1  +  (fi-l,JL  "  f i-l.JL-l5 


(91) 


In  this  case,  the  quantity  w  is  defined  as 

w  »  tan  (y  +  0) 
m 

and  by  geometry 

v  -  (Ar  -  i) /Ax 
P 


The  interpolation  scheme  for  Wp  now  gives  the  expression: 


v 

P 


Wi-1 , JL-1  +  (wi-l,JL  ”  Wi-1,JL-1) 


(92) 


(93) 


(94) 
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Again,  equating  Equation  93  to  Equation  94  and  solving  for  i  gives: 

t  .  (Ar/4x)  -  w1.ltJ^.1 _  (95) 

(  i-l.JL  "  Wi-1,JL-1)/Ar  +  (1/Ax) 

Applying  the  computed  values  of  the  interpolation  length  to  the 
respective  interpolation  equation  (either  Equation  86  or  Equation  91) 
gives  the  proper  value  of  the  desired  flow  variable  at  the  boundary  point 
(i,  JL). 

3.  THE  DOWNSTREAM  BOUNDARY 

The  downstream  boundary  (EF)  is  placed  in  a  region  where  gradients  in 

the  flow  variables  are  expected  to  exist.  A  quadratic  extrapolation  can  be 

2  2 

used  on  this  boundary  that  lets  3f/3x  and  3  f/3x  be  nonzero,  thus  satis¬ 
fying  this  gradient  condition.  Assuming  a  constant  grid  spacing  Ax  near 
this  boundary,  Taylor  series  expansions  can  be  performed  in  the  following 
manner  from  a  point  (IL,  j)  on  the  boundary  to  the  following  three  points 
interior  to  the  boundary: 
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2  2  i  3  i 

.  ax*  rf  , 

IL,j  2  ax2|lL,3  °(  3x3'IL*j> 


*  a.  o.  if  .  -A  2  32f  .  rt,.  3  33f 

f IL, j  2&x  ax  IL, j  2Ax  a„2  IL, j  0(Ax  »_3 


f  it 

£il,j  +  J“x  ax1 


ax 

9 A X2  a2f 


3x 

3 

3  aJfx 


(96) 


XL.j’  (»> 

n.,3  +  ~ i T  ~^i |il,j  +  0<lx  liL.j'  <98> 


If  the  assumption  is  made  that  the  last  term  in  each  equation  can 
be  neglected,  i.e. 


0 

axJ 


(99) 
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then  Equations  96  through  98  can  be  solved  simultaneously  to  give  the 
following  expression  for  the  boundary  point  (IL,  j )  in  terms  of  the 
interior  points: 

f IL, j  “  3  f IL-1 , j  "  3  fIL-2,j  +  f  IL-3, j  (100) 

This  condition  works  well  as  long  as  large  pressure  gradients  do 

not  exist  at  this  boundary,  as  in  the  case  when  a  normal  shock  wave 

3  3 

exits  the  boundary.  If  this  does  occur,  the  term  involving  3  f/3x  is 
no  longer  negligible.  Equation  100  then  can  become  numerically  unstable. 

Numerical  divergence  did  occur  when  the  previous  extrapolation 
condition  was  applied  to  regions  of  subsonic  flow  present  at  this  out¬ 
flow  boundary.  Therefore,  the  following  first  order,  zero  gradient 
condition  was  applied  at  points  (IL,  j)  when  the  Mach  number  at  points 
(IL-1,  j)  was  found  to  be  subsonic: 

f IL,  j  "  f IL-1 ,  j  (101  ) 

4.  THE  CENTERLINE 

The  centerline  boundary  (AF  in  Figure  6)  is  a  line  of  symmetry  with 
no  mass  or  energy  flux  across  it.  The  vertical  velocity  component 
condition  (Equation  36)  is  applied  by  setting: 

Vj  j  ■  O  for  1  <_  i  <_  IL  (102) 

The  symmetry  conditions  for  the  density  and  u  component  of  velocity 
(Equations  37  and  38)  which  are  valid  only  at  the  centerline  are  applied 
in  the  following  manner.  Taylor  series  expansions  are  performed  for  p 
and  u  from  a  centerline  boundary  point  (i,l)  to  points  (i,2)  and  (i,3) 
which  are  distances  Ar  and  KAr,  respectively,  above  the  centerline 


46 


AFWAL-TR-81 -31 61 


boundary.  The  series  expansions  to  these  points  give  the  following 
equations  for  density: 


°i,2  ■  “i.l  +  4t  +  0^3 


»1.3  -  *1.1  +  t 


2  2  2 
R Ar  ac 


.  k  nr  api  n  3  33d  , 

i  1  +  5 - 9  \i  i  +  °(Ar  — r  .  ,) 

2  ar2|i.l  3t3  i»1 


(103a) 


(103b) 


If  the  centerline  symmetry  condition  (Equation  37)  is  applied  and 
the  higher  order  term  in  each  equation  is  neglected,  these  two  equations 
can  be  solved  simultaneously  to  obtain: 


K2  pi,2  -  pl,3 
<K2  -  1) 


(103c) 


A  similar  expression  can  be  obtained  for  the  u  component  of  velocity: 


V2 

K  uj,2  ~  ui,3 
<K2  -  1) 


(103d) 


The  previous  extrapolation  boundary  conditions  for  the  density  and 
the  horizontal  velocity  component  were  applied  only  to  regions  of  super¬ 
sonic  flow  at  the  centerline.  Undesirable  pressure  wiggles  occurred  if 
these  conditions  were  applied  to  regions  of  subsonic  flow.  The  following 
first  order  zero  gradient  condition  was  applied  at  the  centerline  points 
(i,l)  when  the  Mach  number  at  points  (1,2)  was  found  to  be  subsonic 


Pi.l  "  Pi,2 


(103e) 


ui,l  "  ui,2 


( 1 03f ) 


Since  the  v  component  of  velocity  is  zero  on  the  centerline,  this 
boundary  can  be  considered  as  a  streamline.  As  discussed  in  Section  II, 
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the  stagnation  enthalpy  is  then  constant  on  this  boundary.  Since  the 
total  enthalpy  at  the  inflow  boundary  is  known,  then 


H 

o 


1,1 


H  .  .  for  1  <  i  <  IL 

o  1,1  — 


(104) 


Since  the  value  of  u^  ^  has  been  previously  determined  by  Equation  103, 
the  definition  of  the  stagnation  enthalpy  can  be  expanded  to  determine 
the  value  of  temperature  at  each  boundary  point  (i,l): 


li,l 


Ti,i+  '  (“i,i>2>/2cp 


(105) 


The  values  of  the  primitive  variables  p,  u,  v,  and  T  have  now  been 
determined  for  each  centerline  boundary  point,  so  that  the  required 
values  of  the  conservative  variables  can  be  computed  along  this  boundary. 

5.  THE  NOZZLE  WALLS 

The  nozzle  walls  (BG,  GH,  and  CH  in  Figure  6)  are  treated  as 
no-slip,  impermeable  boundaries.  The  no-slip  condition  is  imposed  on 
the  three  wall  faces  by  imposing  the  following  conditions  (see  Figure  17): 

Inner  wall 


u 


i,JWI 


Vi , JWI 


-  0  for  1  <  i  <  IW 


Outer  wall 


Ui,JVIO  =  Vi,JWO 


0  for  1  <  i  <  IW 


Base  (vertical )  wal  1 


uiw,j  "  viw,j 


0  for  JWI  <  j  <  JWO 


(106a) 


(106b) 


(106c) 
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O  JWO+1 

O  JWO 

O 

o 


o 

O  JWI 

O  JWi-i 
IW+l 


Figure  17.  Finite  Difference  Mesh  Near  the  Nozzle  Walls 
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As  discussed  in  Section  II,  a  constant  wall  temperature  is  imposed 
on  the  nozzle  walls.  This  condition  is  applied  simply  as: 

T  =  T  «=  T  for  1  <  i  <  IW  (107a) 

i.JWI  i.JWO  W  - 

and 

T  »  T  for  JWI  <  j  <  JWO  (107b) 

IW ,  j  w 

A  first  order  pressure  gradient  condition  derived  from  Equation  48 
is  applied  on  each  wall.  This  is  imposed  as: 


i ,  JWI 

"  Pi, JWI-1  for  1 

<  i 

<  IW 

i.JWO 

Pi , JWO+1  f0r  1 

<  i 

<  IW 

iw.j  ■ 

PIW+1 , j  for  JWI 

<  j 

<  JWO 

(108a) 

(108b) 

(108c) 


Since  the  points  (IW,  JWI)  and  (IW,  JWO)  are  positioned  at  sharp 
corner  points,  the  simple  pressure  condition  applied  in  Equation  108  is 
not  applicable.  An  averaging  scheme  was  therefore  used  to  allow  the 
pressure  to  adjust  at  the  corners.  This  averaging  is  applied  as: 


P  =  (r*  +  p  i  /2 

IW.JWI  '  IW,  JWI-1  rIW+l,  JWl" 

and 

P  -  (p  +  p  )/2 

IW.JWO  v  IW,  JWO+1  IW+1,  jwo" 


(109a) 

(109b) 


The  primitive  variables  u,  v,  P,  and  T  have  now  been  defined  on 
the  wall  boundaries.  The  required  values  of  the  conservative  variables 
can  then  be  computed  on  this  boundary. 


6.  INITIAL  CONDITIONS 

As  discussed  in  Section  II,  initial  values  of  the  conservative 
variables  must  be  imposed  over  the  computational  domain.  Since  the 
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incoming  flow  variables  are  fixed  in  time,  they  are  initially  imposed 
over  the  complete  computational  domain  as  follows: 


U°  =  U°  ,  j  <  JWI  and  JWO  <  j 

J  J 


(no) 


The  value  of  the  u  component  of  velocity  in  the  wake  region  (JWI  £  j 
<_  JWO)  would  be  zero  from  the  input  profile.  Therefore,  the  u  com¬ 
ponent  in  the  wake  is  set  to  grow  exponentially  using  the  following 
equation : 


Xj 


k„(l  -  e[  (xIW,j  ‘  xi,j)/yBJ) 


(HD 


where 
k0  - 


x  X  *  ~  r! 


JWI-1 


) 


(112) 


for  IW  <  i  ^  IL  and  JWI  <  j  <  JWO.  Use  of  this  scheme  allows  the  velocity 
in  the  far  wake  to  be  close  to  that  of  the  two  streams,  thus  accelerating 
convergence. 

Once  an  initial  case  for  the  nozzle  flowfield  had  numerically 
converged  to  a  valid  solution,  each  succeeding  case  was  initialized  by 
imposing  a  new  jet  input  profile  on  the  inflow  boundary  of  the  preceding 
converged  solution.  This  technique  allowed  the  new  solution  to  converge 
at  a  much  greater  rate,  since  the  subsonic  recirculation  zone  in  the 
near  wake  was  already  in  existence  from  the  previous  solution. 
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SECTION  V 

TURBULENCE  MODELING 

The  experimental  tests  used  as  a  basis  for  the  computational  solutions 

were  conducted  at  a  Reynolds  number  of  2.2  x  10^,  based  on  the  ogive  body 

length  and  the  external  flow  conditions.  The  external  flow  in  the  region 

of  the  nozzle  is  thus  expected  to  be  of  a  fully  turbulent  nature.  Reynolds 

5  6 

numbers  in  the  interior  jet  flow  covered  a  range  from  1  x  10  to  1.7  x  10  , 
based  on  nozzle  exit  conditions  and  the  nozzle  throat  to  exit  plane  length. 
Considering  the  effects  of  compressibi 1 ity  and  the  existence  of  a  favorable 
pressure  gradient  in  the  divergent  portion  of  the  nozzle,  a  transition 

5 

Reynolds  number  of  5  x  10  was  used  to  determine  the  condition  at  which  the 
jet  flow  possessed  turbulent  characteristics  (Reference  27). 

The  turbulent  nature  characteristic  of  these  flows  can  be  accounted 
for  in  the  computational  solution  by  a  variety  of  eddy  viscosity  models 
ranging  from  locally  dependent  algebraic  models  to  the  more  complex  higher- 
order  closure  models  such  as  the  turbulent  kinetic  energy  methods. 

Although  the  higher  order  methods  can  account  for  the  time  history  of 
the  turbulence  in  a  flow,  they  require  that  accurate  initial  profiles  of 
the  turbulent  shear  stress  be  known  or  reliably  calculated  (Reference  28). 
If  this  initial  profile  condition  cannot  be  satisfied,  then  this  type  of 
prediction  method  cannot  be  effectively  utilized.  Since  this  proved  to 
be  the  case  for  the  jet  problem  under  consideration,  locally  dependent 
eddy  viscosity  models  were  carefully  applied  over  the  computational 
domain. 

As  shown  in  Figure  18,  the  computational  domain  contains  three  dis¬ 
tinct  regions  in  which  various  eddy  viscosity  models  are  applied.  These 
regions  consist  of  an  area  containing  boundary  layers,  a  far  wake  region 
downstream  of  the  nozzle  exit,  and  a  near  wake  region  close  to  the  nozzle 
exit  plane. 

1.  BOUNDARY  LAYER  MODEL 

In  the  first  region,  the  dominating  flow  features  are  the  boundary 
layers  along  the  nozzle  walls.  Since  the  experimental  boundary  layer 
thicknesses  are  at  least  an  order  of  magnitude  smaller  than  the  nozzle 
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Figure  18.  Eddy  Viscosity  Model  Domains 


radius,  a  two-dimensional  turbulence  model  was  judged  to  be  sufficient 
for  the  axi symmetric  cases.  The  eddy  viscosity  model  applied  in  this 
region  is  a  two  layer  Cebeci-Smith  model  (Reference  29).  The  inner  layer 
of  this  model  accounts  primarily  for  the  laminar  sublayer  adjacent  to 
the  wall,  with  the  outer  layer  accounting  for  the  remainder  of  the 
boundary  layer  region. 

The  expression  for  the  inner  model  is  based  on  Prandtl's  mixing- 
length  theory,  which  can  be  written  as: 


p 


3r 


(113) 


where  ut  is  the  local  tangential  velocity  parallel  to  the  wall  surface, 
and  rn  is  the  normal  distance  measured  from  the  wall.  The  mixing  lei.gth 
in  this  model  is  adapted  from  Van  Driest's  sublayer  model,  and  is  ex¬ 
pressed  as: 


-  r  fpi 
n  v 


0 . A  r  (1  - 


26p 


014) 
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This  inner  model  assumes  the  flat  ^late  pressure  condition  (dp/dx  =  0), 
but  could  be  modified  to  account  for  a  pressure  gradient  in  the  direction 
parallel  to  the  wall  within  the  sublayer. 

The  outer  region  eddy  viscosity  model  consists  of  a  Clauser-type 
displacement  thickness  model  defined  by  the  equation: 

c0  “  0.0168  pug  6  (115) 

where  ug  is  the  appropriate  tangential  velocity  at  the  boundary-layer 
edge  and 


is  the  incompressible  displacement  thickness.  This  model  also  includes 
Klebanoff's  intermi ttency  factor  defined  by  the  following  equation: 

Yj  -  [1  +  5.5.  (y/6)6]"1  (117) 

The  inner  and  outer  regions  of  each  boundary  layer  are  defined  by 
the  requirement  that  the  eddy  viscosity  remain  continuous  across  the 
entire  layer.  This  is  accomplished  by  applying  the  inner  model  outward 
from  the  wall  until  =  c0  at  a  value  r  .  The  outer  model  is  then  applied 
from  rc  outward  across  the  remainder  of  the  flowfield  in  the  boundary 
layer  region.  Figure  19  shows  a  typical  eddy  viscosity  profile  across 
this  region. 

2.  FAR  WAKE  MODEL 

In  the  region  downstream  of  the  nozzle  exit,  the  initial  boundary 
layers  on  the  nozzle  walls  merge  to  form  a  shear  layer  containing  an 
imbedded  wake  region.  This  region  in  the  flowfield  can  be  further  divided 
into  two  separate  areas:  the  near  wake  region  close  to  the  nozzle  exit 
that  contains  flow  features  such  as  the  corner  expansions,  a  "deadwater" 

zone,  recompression  shocks,  and  the  far  wake  region  further  downstream  > 

where  the  flow  eventually  tends  to  a  similar  free  shear  layer  type  of  flow. 

The  eddy  viscosity  model  for  the  far  wake  region  will  be  discussed  in  this 
section,  and  then  be  extended  for  the  near  wake  region  in  the  next  section. 
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The  eddy  viscosity  in  the  far  wake  region  is  in  the  form  of  the 
following  Prandtl  mixing  length  model: 

e  -  P*2|u| 


(118) 


where  <o  is  defined  as  the  vorticity: 


fiu  fiv 


Sr  fix 


(119) 


and  the  mixing  length  i  uses  the  same  formulation  as  that  of  Dash,  et  al. 
(References  30  and  31)  in  their  wake  mixing  length  model  for  the  core 
region  of  a  coflowing  nozzle: 

l  -  0.065  (120) 

In  this  model,  6w  is  the  representative  thickness  of  the  shear  layer  at 
any  axial  position  in  the  wake.  This  model  accounts  for  the  variation 
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in  eddy  viscosity  across  the  mixing  layer  through  its  dependence  on  the 
local  value  of  vorticity.  As  in  the  eddy  viscosity  model  utilized  by 
Baldwin  and  Lomax  (Reference  32),  the  vorticity  profile  across  the  mixing 
layer  can  be  utilized  to  determine  the  thickness  parameter  6  .  This 
eliminates  the  somewhat  complicated  process  of  finding  the  outer  edges  of 
the  shear  layer  based  on  velocity  profiles  for  each  axial  position  in 
the  computational  flowfield. 

For  the  axi symmetric  coflowing  nozzle  cases  that  were  solved 
numerically,  the  maximum  absolute  value  of  vorticity  in  the  wake  was 
found  to  be  in  the  following  range: 

1  x  105  <  |u  |  <  1  x  106  sec"1  (121 ) 

1  max 1 

The  cutoff  value  used  to  define  a  representative  edge  of  the  mixing 
layer  was : 

I uedge  I  -  1  x  104  sec"1  O22) 

This  value  gave  a  reasonable  value  of  6  as  shown  by  a  typical  vorticity 
profile  in  Figure  20.  The  absolute  value  of  vorticity  typically  dropped 
to  less  than  one  percent  of  [w  x|  within  one  gridpoint  outside  of  the 
cutoff  point. 

A  flat  plate  validation  case  was  computed  using  the  far  wake  model 
to  check  its  accuracy  in  a  known  turbulent  flowfield.  The  data  of 
Toyoda  and  Hiriyama  (Reference  33)  for  a  flat  plate  in  turbulent  flow  at 
a  Mach  number  of  1.6  was  used  as  the  basis  for  a  computational  solution. 
The  velocity  defect  in  the  wake  obtained  both  experimentally  and 
numerically  is  shown  in  Figure  21.  As  shown  by  this  figure,  the  results 
generated  by  the  numerical  turbulence  model  compared  very  well  with  the 
thin  flat  plate  data.  Further  details  of  this  computation  are  listed  in 
Appendix  C. 
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3.  NEAR  WAKE  MODEL 

The  accuracy  of  a  Prandtl  type  mixing  length  turbulence  model  is 
substantially  dependent  on  the  use  of  length  scales  that  are  truly 
representative  of  the  flow  in  a  given  region.  In  the  far  wake  region 
where  a  single  mixing  layer  exists,  the  previous  definition  of  the  length 
scale  involving  6  is  valid.  However,  for  the  case  of  a  nozzle  with 

W 

a  thick  base  annulus,  several  length  scales  need  to  be  defined  in  the 
near  wake  region.  As  shown  in  Figure  22,  the  existence  of  the  subsonic, 
recirculating  "dead  water"  region  adjacent  to  the  nozzle  base  wall 
complicates  the  flow  simulation.  In  the  near  wake,  the  length  scales 
must  transition  from  the  appropriate  boundary  layer  thickness  at  the 
nozzle  exit  plane  to  the  single  mixing  layer  thickness  (6  )  that  exists 
in  the  far  wake  region. 

This  transition  is  accomplished  using  the  following  procedure.  The 
exterior  edges  of  the  mixing  layers  in  the  near  wake  are  determined  using 
the  vorticity  limits  previously  defined  for  the  far  wake.  The  interior 
edges  of  the  dual  shear  layers  are  then  defined  by  the  Mach  0.5  contour 
line  surrounding  the  "dead  water"  region  as  shown  in  Figure  22.  The 
vorticity  turbulence  model  defined  by  Equations  118  and  120  is  then  applied 
in  the  near  wake  region  from  the  end  of  the  boundary  layer  zone  to  the 
start  of  the  initial  mixing  zone  with 

6  =  5 o  in  the  free  stream  flow 

6  =  6 .  in  the  jet  stream  flow 

6  =  0.5  (<5q  +  <5.j)  inside  the  Mach  0.5  contour  of  the  "dead  water" 
region 

An  initial  mixing  zone,  one  base  height  in  length,  is  used  to  smoothly 
adjust  the  thickness  of  the  dual  shear  layers  (6Q  and  6^)  that  exist  in 
the  expansion  and  recompression  zones  of  the  near  wake  region  to  that  of 
the  single  shear  layer  (&)  in  existence  further  downstream.  The  follow- 

W 

ing  exponential  equation  is  applied  in  this  region: 


6(x)  •  6 

w 


-6 

o  or 


(123) 
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The  point  xm  is  centered  in  the  mixing  region.  It  was  determined  that 
the  position  of  the  midpoint  xm  must  be  specified.  If  allowed  to  float, 
the  near  wake  eddy  viscosity  region  stretched  to  unrealistic  lengths 
and  gave  erroneous  values  of  nozzle  base  pressure.  Using  the  very 
slightly  underexpanded  experimental  case  (P./P^  =  1.03)  as  a  basis,  the 

J 

midpoint  xm  was  varied  to  obtain  its  effect  on  nozzle  base  pressure. 

As  shown  in  Figure  23,  the  value 

yn>  -  2-° 

gave  the  best  agreement  with  the  experimental  data.  This  value  of  xm  was 
then  fixed  for  the  entire  series  of  flow  calculations  at  the  various 
pressure  ratios.  The  variation  in  nozzle  base  pressure  shown  in  Figure  23 
illustrates  the  dependence  of  base  pressure  on  the  magnitude  of  the  eddy 
viscosity  applied  in  the  near  wake  region. 

A  two-dimensional  wedge-flat  plate  validation  case  was  computed  in 
order  to  obtain  the  accuracy  of  the  turbulence  model  in  the  near  wake 
and  transition  zone  to  the  far  wake.  The  data  of  Rom,  Seginer,  and 
Kronzon  (References  34  and  35)  for  a  one  centimeter  thick  wedge-flat 
plate  in  turbulent  flow  at  a  Mach  number  of  2.25  was  used  as  the  basis 
for  a  computational  solution. 

All  of  the  near  wake  features  were  accurately  reproduced  by  the 
computational  solution,  and  are  discussed  in  Appendix  D.  Both  the  static 
pressure  axially  along  the  line  of  symmetry  and  the  pitot  pressure  pro¬ 
files  in  the  near  wake  are  in  good  agreement  with  the  experimental  data. 
The  variation  in  the  mixing  length  (s)  used  in  the  turbulence  model  is 
shown  in  Figure  24.  This  figure  exhibits  the  correct  physical  behavior 
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of  the  growth  of  6  as  the  flow  expands  around  the  corner  and  the  subsequent 
decrease  in  6  as  the  flow  recompresses  in  the  near  wake  region.  The 
exponential  growth  in  the  transition  region  is  also  evident.  This  case 
confirms  that  the  eddy  viscosity  model  used  in  the  coflowing  nozzle  is  a 
reasonable  one  that  should  account  for  the  thrublence  effects  in  a  correct 
manner. 
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SECTION  VI 
NUMERICAL  RESULTS 

This  section  will  discuss  the  numerical  results  of  the  computational 
solutions  using  the  algorithm,  boundary  conditions,  and  turbulence 
modeling  detailed  in  previous  chapters.  The  first  portion  of  this  section 
will  discuss  the  experimental  cases  taken  as  the  basis  for  comparison 
with  the  numerical  solutions  for  the  coflowing  nozzle.  The  next  portion 
discusses  details  involved  in  the  actual  computational  procedure.  The 
last  part  covers  the  comparison  between  the  experimental  and  computational 
solutions,  including  some  analyses  of  the  accuracy  of  the  simulations 
and  discrepancies  between  the  numerical  solutions  and  the  experimental 
data . 

1.  EXPERIMENTAL  DATA  BASE 

As  outlined  in  Section  I,  the  data  of  Bromm  and  O'Donnell  (Reference 
16)  is  used  as  the  experimental  basis  for  this  research  effort.  Super¬ 
sonic  fields  of  flow  generated  experimentally  contain  both  highly  viscous 
flow  regions  as  well  as  shock  structures  ranging  from  weak  regularly 
reflected  shock  waves  to  the  strong  Mach  disc  shock  formation.  Five 
different  experimental  nozzle  pressure  ratio  conditions  are  used  as  the 
basis  for  the  computational  solutions.  Nozzle  base  pressure  measurements 
and  schlieren  photographs  are  the  basis  for  experimental  versus  compu¬ 
tational  comparisons. 

a.  Model 

The  model,  a  stainless  steel  body  of  revolution,  consisted  of 
a  one-inch  diameter  cylindrical  body  with  a  16.25  inch  radius  ogive  nose 
(see  Figure  15).  The  total  length  of  the  model  was  7.5  inches.  The 
model  was  supported  by  a  10%  thick  hollow  side  strut  that  acted  as  a 
conduit  for  the  air  flow  to  the  jet  as  shown  in  Figure  5.  The  effects  of 
this  strut  were  found  to  be  negligible  on  the  flowfield  downstream. 

This  model  was  fitted  with  a  nozzle  which  gave  essentially  isentropic 
flow  with  an  exit  Mach  number  of  3.00.  The  inner  diameter  of  this  nozzle 
at  the  exit  plane  was  0.742  inches,  and  the  length  from  the  nozzle  throat 
to  the  exit  plane  was  1.20  inches.  Four  base  pressure  orifices  were  used 
to  obtain  the  base  pressure  measurements  as  shown  on  Figure  5. 
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b.  External  Flow 

The  external  flow  conditions  were  generated  in  the  NASA  Langley 
9-inch  supersonic  wind  tunnel.  The  free  stream  Mach  number  was  set  at 
1.94,  and  the  free  stream  Reynolds  number  was  fixed  at  2.2  x  10^  based 
on  the  body  length  of  the  model.  A  turbulent  boundary  layer  on  the  model 
at  these  conditions  was  insured  through  the  use  of  a  transition  strip 
near  the  nose.  The  tests  were  conducted  at  a  tunnel  stagnation  pressure 
of  one  atmosphere  (assumed  to  be  2116  psf).  Using  these  given  conditions, 
the  stagnation  temperature  of  the  free  stream  was  calculated  to  be  equal 
to  530. 5°R. 

c.  Jet  Flow 

The  Flow  just  upstream  of  the  jet  exit  plane  was  given  to  be  at 
a  Mach  3,  zero  divergence  angle  condition.  The  total  pressure  in  the  jet 
flow  was  varied  to  obtain  the  desired  nozzle  static  pressure  ratio.  The 
jet  static  pressure  at  the  nozzle  exit  plane  was  not  measured  directly, 
but  was  calculated  using  the  given  nozzle  area  ratio  and  the  jet  total 
pressure.  A  total  temperature  in  the  jet  flow  was  not  given  experimentally, 
but  was  assumed  to  be  equal  to  the  freestream  stagnation  temperature 
(580. 5°R). 

2.  COMPUTATIONAL  DETAILS 

Solutions  were  computed  for  the  axisymmetric  nozzle  at  the  following 
five  nozzle  static  pressure  ratios:  Pj/Pra  =  0-150,  0.251,  0.527,  1.03 
and  1.59.  These  solutions  were  all  performed  on  a  CDC  Cyber  175  digital 
computer  located  at  Wright-Patterson  AFB,  Ohio.  The  average  rate  of  data 
processing  was  0.0015  second  per  grid  point  per  iterative  time  increment. 

In  this  section  further  details  of  these  computations  will  be  discussed. 

a.  Grid  Parameters 

Important  parameters  of  the  computational  grid  including  the 
number  of  grid  points  and  axis  length  utilized  in  both  the  axial  and 
radial  directions  are  listed  in  Table  1  for  each  nozzle  pressure  ratio 
condition.  The  value  of  the  axial  field  length  includes  a  length 
increment  of  0.4  upstream  of  the  nozzle  exit  plane,  with  the  exit  plane 
at  a  value  of  x/rjet  =  0.0.  A  compact  45  x  45  point  grid  was  used  in  the 
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TABLE  1 


COMPUTATIONAL  GRID  PARAMETERS 


IL  JL 


£L 


0.150 

45 

45 

3.4 

3.0 

0.251 

45 

45 

3.4 

3.0 

0.527 

57 

45 

4.4 

3.0 

1.03 

57 

45 

6.4 

3.0 

1.59 

57 

45 

8.4 

4.0 

two  most  highly  overexpanded  cases  (Pj/P^  =  0.150,  0.251).  Twelve 
additional  grid  points  were  then  added  in  the  axial  direction  downstream 
of  the  original  grid  to  form  a  57  x  45  grid  for  the  next  case  where 
Pj/P<x>  =  0.527  (Figure  25).  This  methodology  gave  a  consistent  cell 
length  in  the  axial  direction  for  each  of  these  three  cases  in  which 
both  a  regular  shock  reflection  (P./P^  =  0.527)  and  Mach  discs  (P./P^  = 
0.251,  0.150)  should  occur.  The  computational  domain  was  then  stretched 
using  the  57  x  45  point  grid  to  the  degree  necessary  to  cover  the 
phenomena  of  interest  for  the  two  underexpanded  cases. 

Minimum  grid  spacing  in  both  the  axial  and  radial  directions 
occurred  adjacent  to  the  nozzle  walls,  and  was  set  at: 

Ax/rjet  =  Ar/rjet  =  °-030 

The  patched  exponential  stretching  outlined  in  Section  III  was  then 
applied  to  form  each  grid  initially  for  the  various  cases.  The  adaptive 
grid  procedure  was  applied  in  the  radial  direction  during  the  solution 
procedure  as  discussed  in  Section  III  and  Appendix  B  to  obtain  the  final 
grid  geometry  for  each  case. 

b.  Coarse  Grid  Effects  on  Boundary  Layer  Resolution 

The  numerical  solution  of  the  Navier-Stokes  equations  can  involve 
significant  truncation  errors  in  regions  containing  high  velocity  gradients 
such  as  within  turbulent  boundary  layers  when  relatively  coarse  compu¬ 
tational  grids  are  employed.  Errors  in  computed  velocity  gradients 
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x/r 


jet 


Figure  25.  Adaptive  Finite  Difference  Mesh,  P./P^  =  0.527 

vj 

involved  in  shear  force  terms  at  the  wall  can  result  in  erroneous  values 
of  the  pressure  gradient  along  the  wall.  The  extent  of  these  errors  can 
be  realized  by  assuming  the  existence  of  a  turbulent  boundary  layer 
possessing  a  velocity  profile  in  the  following  form  (Reference  36): 


+ 

u 


+ 

u 


y  i 


[.50  ln(y+)  +  5.10  11  <  y+  <  y+BL  e 


(125a) 

(125b) 


where 

u+  -  u^?  .w  and  y+  -  ^  P 

Since  the  dominant  term  in  the  shear  stress  is  the  gradient 
a  comparison  of  the  value  obtained  using  this  profile  can  be  made 
wi  .h  that  obtained  using  the  numerical  algorithm.  MacCormack's  algorithm 
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computes  this  gradient  using  a  first  order  finite  difference  in  the 
direction  of  the  sweep  for  the  predictor  step.  At  the  wall  this  gradient 
is  computed  as: 


Au 

Ay 


w 


Ui,j^l/Ayv 


(126) 


where  jw+1  is  the  first  grid  point  above  the  wall,  and  Ay  is  the  grid 

w  0 

spacing  adjacent  to  the  wall.  If  a  nominal  Reynolds  number  of  1x10 
based  on  body  length  is  applied  to  the  previous  turbulent  velocity 
profile,  the  ratio  of  actual  wall  velocity  gradient  given  by  the  previous 
profile  to  that  computed  using  Equation  126  can  be  displayed  as  a  function 
of  the  ratio  of  u  velocity  component  at  the  first  point  above  the  wall  to 
the  freestream  velocity  as  shown  in  Figure  26.  For  grid  spacing  sub¬ 
stantially  greater  than  the  sublayer  ( ui /ue  =  0.475),  excessive  error 
in  the  computed  velocity  gradient  occurs.  Since  in  the  nozzle  solutions 
the  previously  stated  minimum  grid  spacing  gives  values  of  (u-j/u^l^  =  0.72 
and  (u]/ue)ljet  =  0.90,  a  correction  is  needed.  This  is  accomplished  by 
applying  the  value  of  Cf  given  by  a  boundary  layer  analysis  at  the  inflow 
boundary  along  the  nozzle  wall  when  j  =  jw,  instead  of  using  the  standard 
differencing  procedure  in  the  algorithm.  This  procedure  has  resulted  in 
smooth  pressure  profiles  near  the  inflow  boundary  instead  of  slight 
pressure  jumps  previously  observed  in  coarse  mesh  cases.  This  concept  is 
analogous  to  the  wall  functions  described  by  Launder  and  Spaulding 
(Reference  37)  and  employed  by  Peery  and  Forester  (Reference  38). 


c.  External  Flow  Parameters 

Values  for  the  external  flow  variables  were  input  numerically  by 
aDplying  the  experimental  values  of  M^,  P0oo,  T0oo,  and  the  boundary  layer 
profiles  generated  for  use  on  the  inflow  boundary.  A  value  for  the 
turbulent  skin  friction  coefficient  of  =  0.00286  was  obtained  from 
the  boundary  layer  input  profiles  and  was  applied  on  the  external  nozzle 
wall  as  discussed  in  the  last  section.  A  wall  temperature  of  551 °R  was 
calculated  using  the  analysis  in  Appendix  A.  The  wall  temperature  was 
given  a  constant  >'alue  for  all  five  cases  since  the  analysis  showed  that 
this  temperature  should  not  vary  more  than  one  degree  over  the  range  of 
simulated  flow  conditions. 
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d.  Jet  Flow  Parameters 

The  numerical  jet  field  of  flow  was  determined  by  applying  the 
experimental  values  of  the  jet  Mach  number  and  total  temperature  given  in 
the  previous  section  along  with  the  values  of  the  jet  total  pressure  and 
skin  friction  coefficient  listed  in  Table  2.  Since  the  nozzle  wall 
temperature  was  considered  to  be  constant,  the  value  for  the  internal 
nozzle  wall  was  also  set  at  551°R.  As  shown  in  the  table,  the  state  of  the 
jet  boundary  layer  in  all  of  the  cases  except  the  lowest  nozzle  pressure 
ratio  condition  (Pj/P^  =  0.150)  was  considered  to  be  turbulent  at  the 
nozzle  exit  plane.  Since  the  case  where  P./P^  =  0.251  possessed  a  Reynolds 
number  below  the  previously  set  transition  point,  a  solution  using  a 
laminar  jet  boundary  layer  was  also  obtained.  At  this  pressure  condition, 
both  the  laminar  and  turbulent  jet  boundary  layers  gave  nearly  identical 
values  for  the  nozzle  base  pressure  and  normal  shock  position.  However, 
the  laminar  case  exhibited  a  mild  shear  layer  oscillation  for  the  duration 
of  the  solution  procedure.  For  this  reason  the  results  of  the  turbulent 
case  are  presented.  It  is  interesting  to  note  that  no  mixing  layer 
oscillations  were  observed  in  the  other  laminar  jet  condition  (P-/P  = 

0.150). 


TABLE  2 

COMPUTATIONAL  JET  FLOW  PARAMETERS 

P./P^,  P  (psf)  Re  xlO  ^  Bound,  Layer  Cf  x  10^ 
_ _ jet _ _ _ Character _ j _ 


0.150 

1636. 

1.61 

Laminar 

1.60 

0.251 

2738. 

2.70 

Turbulent 

3.39 

0.527 

5747. 

5.67 

Turbulent 

2.86 

1.030 

11230. 

11.08 

Turbulent 

2.50 

1.590  . 

17340. 

17.11 

Turbulent 

2.30 

e.  Boundary  Condition  Application 

The  specific  numerical  boundary  conditions  given  in  Section  IV 
were  applied  to  the  computational  domain  in  order  to  achieve  solutions 
for  the  axisymmetric  nozzle.  No  significant  problems  were  encountered 
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with  the  application  of  the  upstream  boundaries,  the  upper  boundary,  or 
the  nozzle  wall  boundaries.  Difficulties  were  encountered  with  both  the 
downstream  boundary  and  the  centerline  boundary  when  a  significant  portion 
of  the  flow  along  each  boundary  was  subsonic  in  nature.  This  condition 
only  occurred  for  the  highly  overexpanded  case  (Pj/Poj  =  0.150),  where 
a  substantial  embedded  subsonic  region  exists  downstream  of  the  Mach 
disc  shock  structure. 

The  quadratic  extrapolation  given  by  Equation  100  produced  very 
reasonable  results  for  those  cases  where  the  outflow  was  either  entirely 
supersonic,  or  only  subsonic  at  a  very  few  points  in  the  nozzle  wake. 

For  the  case  where  a  substantial  area  of  subsonic  flow  existed  behind  a 
Mach  disc  structure  and  extended  to  the  downstream  boundary,  numerical 
divergence  occurred  when  Equation  100  was  applied.  A  second  order  zero 
gradient  condition  was  then  applied  in  regions  of  subsonic  flow  along 
this  boundary.  Application  of  this  condition  did  not  produce  numerical 
divergence,  but  did  give  unrealistic  pressure  jumps  at  this  boundary. 

A  first  order  zero  gradient  condition  given  in  Equation  101  was  then 
successfully  applied  to  subsonic  regions  on  this  boundary  with  reasonable 
results. 

An  almost  identical  situation  occurred  along  the  centerline 
boundary  for  subsonic  regions  containing  fairly  strong  radial  flow 
gradients  close  to  the  centerline.  Both  the  extrapolation  condition  given 
by  Equations  103c,  d  and  a  second  order  zero  gradient  condition 
produced  unrealistic  radial  oscillations  in  the  numerical  solution 
(called  wiggles)  within  these  regions  of  subsonic  flow.  The  application 
of  a  first  order  zero  gradient  condition  given  by  Equations  103e,  f  helped 
reduce  these  oscillations  to  achieve  a  reasonable  solution. 

Two  solutions  were  also  computed  where  first  the  downstream 
boundary,  and  subsequently,  the  upper  boundary  were  repositioned  greater 
distances  from  the  nozzle  as  discussed  in  Appendix  E.  No  changes  were 
detected  in  either  the  shock  structure  or  the  nozzle  base  pressure 
coefficient.  These  test  cases  validate  the  effectiveness  of  these 
boundary  conditions  at  the  positions  utilized  in  the  actual  nozzle 
solutions. 
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f.  Convergence 

As  stated  in  Section  III,  the  numerical  solutions  were  either 
initially  started  using  only  the  boundary  layer  profiles  across  the 
computational  domain,  or  restarted  from  a  previous  solution  by  applying 
new  input  profiles  at  the  jet  inflow  boundary.  Solution  times  based 
upon  the  convergence  criteria  discussed  in  Section  III  varied  significantly 
for  the  two  methods  of  initial  startup.  Solution  times  on  the  Cyber  175 
using  only  the  boundary  layer  profiles  to  form  the  initial  conditions 
were  approximately  3.0  hours  for  the  45  x  45  point  mesh  and  3.8  hours  for 
the  57  x  45  point  mesh.  Solution  times  for  cases  restarted  from  previous 
solutions  were  approximately  1.7  hours  for  the  45  x  45  mesh  and  2.1  hours 
for  the  57  x  45  mesh.  The  large  difference  in  these  solution  times  is 
mainly  attributed  to  the  length  of  time  required  for  the  subsonic  re¬ 
circulation  region  in  the  near  wake  to  form  and  achieve  a  steady  state 
condition. 

3.  COMPARISON  WITH  EXPERIMENTAL  DATA 

Comparisons  between  the  numerical  solutions  and  the  experimental 
data  can  be  made  both  qualitatively  and  quantitatively.  Figures  27 
through  31  give  a  good  visual  comparison  between  the  numerical  solutions 
depicted  as  Mach  number  contours  and  the  experimental  schlieren  photo¬ 
graphs.  In  these  figures  the  computed  solutions  above  the  centerline 
were  reflected  to  give  a  total  nozzle  flowfield  to  compare  with  the 
schlieren  photographs.  All  features  typical  of  afterbody  types  of  flows 
such  as  the  shock  structure  internal  to  the  jet  core  flow,  external 
recompression  shocks,  and  shear  layer  development  are  readily  discernible 
and  in  very  good  agreement  with  the  experimental  data. 

As  shown  by  the  previous  figures,  the  pressure  condition  at  which 
the  jet  flow  shock  structure  transitions  from  a  regular  reflection  on 
the  centerline  to  the  Mach  disc  formation  lies  between  the  two  cases 
where  P .  / P^  =  0.527  and  P./Pot  =  0.251.  Although  the  shock  structure 

J  «J 

near  the  centerline  appears  similar  in  the  computational  solutions  for 
these  two  cases,  an  enlargement  of  this  region  as  shown  in  Figure  32 
reveals  several  differences.  The  shock  strength  (related  to  the  Mach 
number  jump  across  the  shock)  is  much  greater  in  both  of  the  strong  Mach 
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disc  cases  (P./P^  =  0.150  and  0.251)  then  in  the  regularly  reflected  case 
J 

(P./P^  =  0.527).  The  sonic  lines  in  this  region  are  displayed  as  dashed 

J 

lines  in  Figure  32  in  order  to  easily  identify  regions  of  subsonic  flow. 
Both  Mach  disc  cases  contain  areas  of  subsonic  flow  downstream  of  the 
shock  along  the  centerline,  whereas  the  minimum  Mach  number  behind  the 
regularly  reflected  shock  is  approximately  equal  to  1.65. 

A  check  was  made  on  the  solution  to  the  highly  overexpanded  case 
(P./Pc  =  0.150)  to  determine  if  the  numerical  solution  correctly  simulated 

J 

the  flow  conditions  across  the  strong  normal  shock  in  the  region  near  the 
centerline.  Since  the  v  velocity  components  are  very  small  near  the 
centerline,  a  one-dimensional  analysis  based  on  the  Rankine-Hugoniot 
relations  across  normal  shocks  can  be  applied.  As  shown  in  Table  3,  the 
computational  solution  was  within  2%  of  the  exact  one-dimensional  analysis. 


TABLE  3 


COMPARISON  BETWEEN  A  1-D  ANALYSIS  AND  THE  COMPUTATIONAL 
SOLUTION  ACROSS  THE  MACH  DISC  FOR  P^/P^  =  0.150 


M. 


M„ 


VIl 


h!h. 


VTi 


Exact  (1-D) 

3.00 

.475 

10.33 

3.857 

2.679 

Computational 

3.00 

.442 

10.52 

3.906 

2.695 

Z  Error 

— 

1.1 

1.8 

1.3 

0.6 

Several  other  phenomena  associated  with  afterbody  flows  are  evident 
in  Figure  33,  which  displays  computed  velocity  profiles  at  given  axial 
stations  for  the  large  Mach  disc  case.  The  separated  "deadwater  zone" 
of  recirculating  flow  is  readily  apparent  in  the  near  wake  region,  as  is 
the  development  of  the  near  wake  to  a  far  wake  velocity  profile.  The 
existence  of  the  strong  Mach  disc  near  the  centerline  is  very  evident, 
and  the  flow  in  the  subsonic  core  region  behind  the  shock  accelerates  in 
the  correct  manner  to  a  slightly  supersonic  condition  at  the  outflow 
boundary. 

A  closer  look  at  the  near  wake  region  is  shown  in  Figure  34  for  two 
of  the  computational  cases.  This  figure  illustrates  the  change  in  the 
shape  of  the  "deadwater  region"  from  a  predominantly  symmetric  nature 
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Figure  33.  Computed  Velocity  Profiles,  P./P  =  0.150 
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at  P./P,,,,  =  1.03  to  one  with  an  asymmetric  nature  at  P./P^  =  0.150.  In  this 
figure  the  dashed  lines  denote  the  dividing  streamline  and  the  streamlines 
through  the  stagnation  point  in  the  near  wake  flow  for  each  case.  As 
shown  in  Figure  35,  the  dividing  streamline  moves  toward  the  inner  wall 
of  the  nozzle  as  the  total  pressure  in  the  jet  is  decreased.  Although 
the  stagnation  point  in  the  near  wake  region  moves  radially  as  the  jet 
stagnation  pressure  is  changed,  it  remains  in  a  relatively  constant 
position  axially  for  the  cases  computed. 

Quantitative  comparisons  between  the  experimental  data  and  the  com¬ 
puted  solutions  are  based  primarily  on  two  parameters:  the  axial  distance 
along  the  centerline  from  the  nozzle  exit  plane  to  the  point  of  reflection 
of  the  incident  shock  wave  at  the  line  of  symmetry,  and  the  \  ,lue  of  the 
nozzle  base  pressure  coefficient.  This  reflection  length,  along  with  the 
type  of  shock  reflection  (either  strong  or  weak),  is  a  good  indication 
that  the  inviscid  flow  features  in  the  jet  core  caused  by  viscous-inviscid 


Figure  35.  Position  of  the  Dividing  Streamline  in  the  Computational 
Nozzle  Solutions 
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interaction  are  properly  simulated.  These  computed  shock  reflection 
lengths  are  obtained  by  examining  the  axial  variation  in  Mach  number 
along  the  centerline  as  shown  in  Figure  36  for  a  typical  case.  As  shown 
in  this  figure,  the  shock  reflection  is  diffused  over  three  cell  lengths, 
with  the  computational  value  of  the  reflection  length  taken  as  being  at 
the  midpoint  of  these  three  cells.  Comparisons  between  the  experimental 
and  computational  values  of  these  reflection  lengths  are  shown  in  Figure  37 
and  Table  4.  Excellent  agreement  was  obtained,  with  the  computational 
results  being  within  2%  of  the  experimental  data. 

An  additional  quantitative  comparison  was  also  made  of  the  computed 
and  observed  Mach  disc  radii  for  the  two  cases  at  which  the  Mach  disc  was 
observed.  This  comparison  is  listed  in  Table  5.  The  computed  Mach  radius 
was  taken  as  the  radial  height  of  the  sonic  line  immediately  behind  the 
strong  shock  wave.  Very  good  agreement  was  obtained  for  the  highly 


Figure  36.  Axial  Variation  in  Computed  Centerline  Mach  Number, 

P./P  =  0.150 

J  00 
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TABLE  4 


COMPARISON  OF  SHOCK  REFLECTION  LENGTHS 


Experimental 


Computational 

xr/r1et(±°*04) 


Z  Error 


0.150 

1.19 

0.251 

1.88 

0.527 

2.98 

1.030 

4.37 

1.590 

5.64 

1.17 

-0.3 

1.91 

+0.5 

2.91 

-1.2 

4.26 

-1.9 

5.53 

-2.0 

overexpanded  case  (Pj/P^  =  0.150)  in  which  a  fairly  large  Mach  disc  occurs. 
Poorer  agreement  was  obtained  for  the  case  that  was  much  nearer  to 
transition  to  a  regular  shock  reflection,  with  a  very  small  Mach  disc 
radius  ( rm/rjet  =  0.17).  The  first  cell  height  adjacent  to  the  centerline 
in  the  computational  solution  possessed  a  value  of  r/r.^  =  0.06,  so  that 
numerical  truncation  error  played  a  large  part  in  the  discrepancy  between 
the  experimental  and  computed  values  of  the  Mach  disc  radial  height  at  this 
pressure  ratio  condition. 

Comparisons  between  the  experimental  and  computational  values  of 
nozzle  base  pressure  coefficient  are  shown  in  Figure  38  and  Table  6. 

Since  the  experimental  data  points  for  the  base  pressure  coefficients 
were  not  obtained  at  the  same  pressure  ratio  values  as  the  Schlieren  data, 
experimental  valuer  for  the  base  pressure  coefficients  in  Table  6  were 
interpolated  from  the  available  data  points  at  the  five  given  nozzle 
pressure  ratios.  Values  of  the  computed  nozzle  base  pressure  were  in 
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Figure  38. 

Base  Pressure  Coefficient  vs  Nozzle  Pressure 

Ratio 

TABLE 

6 

COMPARISON  OF  BASE  PRESSURE  COEFFICIENTS 

V*. 

Experimental* 

Computational 

Z  Error 

p„ (+0.003) 

p, (+0.002) 

a  o 

0.150 

-0.240 

-0.310 

-25.0 

0.251 

-0.280 

-0.299 

-  6.8 

0.527 

-0.265 

-0.276 

-  3.9 

1.030 

-0.237 

-0.245 

-  2.9 

1.590 

-0.226 

-0.219 

+  2.5 
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good  agreement  with  the  experimental  data  (3-7%  error),  with  the  exception 
of  the  highly  overexpanded  case  at  which  Pj/P^  =  0.150.  Figure  38  shows 
that  as  the  pressure  ratio  of  the  nozzle  is  lowered,  the  decreasing  trend 
in  nozzle  base  pressure  reverses  at  a  value  of  approximately  P^/P^  ~  0.18 
and  sharply  increases  as  the  pressure  ratio  is  further  reduced.  This 
sudden  reversal  in  behavior  is  apparently  due  to  flow  separation  in  the 
divergent  portion  of  the  nozzle  which  prevents  the  jet  flow  from  expanding 
fully  to  its  assumed  Mach  3.0  state.  For  pressure  ratio  values  less  than 
Pj/Poo  =  0.135,  the  shock  structure  in  the  jet  core  transitions  back  to  a 
regular  shock  reflection  as  shown  in  Figure  39.  The  schlieren  photographs 

in  Figure  39  indicate  that  for  values  less  than  P-/P  =  0.135,  either  the 

J  00 

jet  Mach  number  is  less  than  1.48  and  thus  cannot  support  a  Mach  disc 
structure  (Reference  8),  or  the  increase  in  nozzle  static  pressure  due  to 
the  reduced  expansion  cannot  produce  the  deflection  angle  in  the  jet  flow 
needed  for  a  Mach  disc  to  occur.  A  non-separated  condition  was  assumed 
by  the  experimental  investigators,  since  their  value  of  P.  was  determined 

J 

using  the  jet  total  pressure  and  the  final  area  ratio  of  the  nozzle. 
Likewise,  the  computational  solutions  assumed  non-separated  Mach  3  flow 
just  upstream  of  the  nozzle  exit  plane.  If  some  separation  did  occur 
and  the  nozzle  flow  did  not  fully  expand  to  a  Mach  3  condition,  a  sub¬ 
stantial  difference  in  base  pressure  could  result. 

This  hypothesis  of  separation  in  the  nozzle  was  partially  confirmed 
by  computationally  solving  a  case  where  the  jet  total  pressure  corresponded 
to  the  attached  case  (P./P^  =  0.150),  but  with  jet  flow  to  a  Mach  number 
of  only  2.60  as  detailed  in  Appendix  E.  A  correct  strong  shock  structure 
was  obtained  computational ly,  and  the  value  of  the  base  pressure  coef¬ 
ficient  increased  to  a  value  of  -0.265.  This  was  in  much  better  agree¬ 
ment  with  the  experimental  data  at  this  condition  (see  Figure  38).  A 
more  accurate  simulation  of  this  pressure  ratio  condition  would  require 
extending  the  computational  mesh  back  to  the  nozzle  throat.  Separation 
in  the  nozzle  could  then  occur  in  a  direct  manner  in  the  numerical 
solution.  Since  this  would  require  extensive  grid  revisions  as  well  as 
additional  computer  resources,  it  was  considered  to  be  beyond  the  scope 
of  this  investigation. 
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Pj  /Pod  =  0.135  Pj/  Pgd  =  0.150 


Figure  39.  Schlieren  Photographs  (Reference  16)  Showing  the  Eventual 
Deterioration  of  the  Mach  Disc  with  Decreasing  Nozzle 
Pressure  Ratio 
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SECTION  VII 

CONCLUSIONS  AND  RECOMMENDATIONS 

A  numerical  method  of  obtaining  solutions  to  the  Navier-Stokes 
equations  for  axisymmetric  nozzles  in  supersonic  external  streams  has  been 
developed  from  a  selection  of  appropriate  techniques.  Based  on  the 
numerical  analysis  and  computational  results  obtained  through  this  study, 
the  following  conclusions  are  drawn: 

1.  The  numerical  solution  of  the  Navier-Stokes  equations  applied  to 
axisymmetric  nozzles  successfully  reproduced  all  of  the  essential 
flow  features  including  boundary  layers,  corner  expansions,  recompres¬ 
sion  snocks,  the  recirculation  region  adjacent  to  the  nozzle  base  wall 
and  the  evolution  of  the  near  wake  to  a  flow  with  far  wake  behavior. 

2.  The  numerical  method  achieved  a  correct  transition  from  regularly 
reflected  shock  waves  at  the  line  of  symmetry  in  the  jet  core  flow  to 
a  strong  Mach  disc  formation  at  the  appropriate  static  pressure  ratio 
condition  of  the  nozzle.  The  subsonic  embedded  region  immediately 
behind  the  Mach  disc  formation  was  simulated  in  a  correct  manner. 

3.  The  application  of  an  adaptive  grid  scheme  in  the  wake  region  of  the 
nozzle  annulus  successfully  positioned  the  fine  mesh  region  of  the 
computational  grid  in  the  wake  region  which  normally  contains  severe 
flow  gradients.  This  allowed  the  accurate  simulation  of  this  high 
flow  gradient  region  while  conserving  numerical  resources. 

4.  The  nozzle  base  pressure  was  heavily  dependent  on  the  eddy  viscosity 
model  applied  in  the  region  of  the  near  wake.  Once  the  model  was 
tuned  for  the  neutrally  expanded  case  (P./P^  =  1.03),  good  agreement 
was  obtained  computationally  for  all  cases  where  the  flow  obeyed  the 
assumption  of  remaining  attached  in  the  divergent  portion  of  the 
nozzle. 

5.  Boundary  conditions  must  be  carefully  formulated  and  applied  in  order 
to  prevent  physically  unrealistic  results  or  numerical  divergence  of 
the  solution.  Both  the  centerline  and  downstream  boundaries  were 
sensitive  where  regions  of  subsonic  flow  occurred  over  a  substantial 
portion  of  the  boundary.  Both  the  quadratic  extrapolation  used  in 
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regions  of  supersonic  flow  and  a  second  order,  zero  gradient  condition 
caused  either  divergence  or  unrealistic  conditions  at  the  boundary 
when  applied  to  regions  of  subsonic  flow.  A  first  order  zero  gradient 
condition  was  used  successfully  in  these  regions  of  subsonic  flow  and 
found  to  be  superior.  Generalizations  about  the  success  of  this  first 
order  boundary  condition  cannot  be  made,  since  the  degree  of  success 
achieved  is  dependent  on  the  specific  numerical  algorithm  applied  to 
the  problem. 

6.  The  final  steady  state  solutions  were  found  to  be  insensitive  to  the 
initial  conditions  applied  over  the  computational  domain.  However, 
the  time  to  converge  to  the  final  solution  was  highly  dependent  on 
the  application  of  specific  initial  conditions.  In  particular,  the 
region  of  subsonic  reci rculation  in  the  near  wake  was  the  last  region 
in  the  solution  domain  to  converge.  Solutions  started  using  only  the 
boundary  layer  profiles  across  the  domain  required  three  hours  to 
converge  on  a  Cyber  175  (45  x  45  point  mesh),  whereas  cases  started 
from  a  previous  different  jet  total  pressure  condition  but  with  an 
established  near  wake  structure  required  only  1.7  hours  to  converge 
(45  x  45  mesh) . 

7.  To  the  author's  knowledge,  this  is  the  first  full  Navier-Stokes 
solution  that  has  accurately  simulated  the  viscous-inviscid  inter¬ 
actions  present  in  a  dual  stream  nozzle  flowfield  at  off-design 
conditions  where  the  strong  Mach  disc  shock  structure  is  present. 
Mikhail  (Reference  8)  previously  was  unsuccessful  in  reproducing  the 
Mach  disc  reflection  in  a  full  Navier-Stokes  solution  due  to  the 
probable  improper  placement  of  the  jet  inflow  boundary  condition 
which  did  not  allow  the  jet  plume  to  expand  to  the  degree  necessary 
to  generate  a  Mach  disc  reflection. 

Based  on  the  numerical  analysis  and  results  obtained  through  this 

study,  the  following  recommendations  are  made: 

1.  The  present  scalar  computer  code  developed  during  the  course  of  this 
investigation  should  be  vectorized  for  use  on  the  new  generation  of 
"supercomputers"  such  as  the  CRAY-1  or  the  Cyber  203.  Although 
present  solution  times  are  on  the  order  of  two  to  four  hours  when 
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run  on  a  Cyber  175  computer,  a  fully  vectorized  version  of  the  present 
computer  code  can  be  expected  to  converge  within  five  minutes  on  a 
CRAY-1  computer  (Reference  39).  This  will  allow  computation  of  more 
complex  nozzle  geometries  and  better  resolution  in  the  boundary  layers 
through  the  application  of  finer  mesh,  while  holding  costs  to  a 
reasonable  level . 

2.  The  present  numerical  solver  should  be  modified  to  include  the  effects 
of  a  calorically  imperfect  gas  with  variable  specific  heat  and  thermal 
conductivity.  This  modification  would  allow  accurate  simulation  of 
the  temperature  dependent  effects  for  hot  exhaust  nozzles  with  gas 
temperatures  less  than  5000°R  (Reference  19).  Only  minor  revisions  to 
the  existing  computer  code  would  be  required  in  order  to  include  these 
temperature  effects. 

3.  After  the  implementation  of  the  two  previous  recommendations,  it  would 
be  desirable  to  incorporate  the  effects  of  species  mixing  into  the 
numerical  solver.  Many  practical  cases  of  interest  -involve  a  jet 
exhaust  gas  with  a  different  species  than  that  of  the  external 
stream.  This  modification  would  require  a  significant  code  revision, 
since  the  addition  of  the  equation  of  mass  diffusion  would  be  required, 
as  would  the  correct  modeling  of  appropriate  mass  diffusion  coef¬ 
ficients  . 
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APPENDIX  A 

NOZZLE  WALL  TEMPERATURE  CALCULATION 

The  nozzle  wall  boundary  condition  applied  during  the  numerical 
solution  procedure  assumes  a  constant  nozzle  wall  temperature  within  the 
computational  domain.  The  relatively  high  thermal  conductivity  of  the 
stainless  steel  nozzle  makes  this  assumption  valid.  This  wall  temperature 
can  be  calculated  by  applying  a  heat  flux  balance  across  the  freestream 
boundary  layer,  the  nozzle  wall,  and  the  jet  boundary  layer  as  shown  in 
Figure  40.  Conduction  of  heat  in  the  axial  direction  is  neglected  due 
to  the  low  temperature  gradients  in  this  direction. 

Since  both  the  freestream  and  the  jet  flow  are  of  a  high-velocity 
nature,  the  unit  heat  flux  for  either  stream  can  be  written  as 
(Reference  40): 

*  hi  (TaWl  “  V  (A-l ) 


where  h^  is  the  heat  conductance  of  each  flowstream.  The  adiabatic  wall 
temperature,  T  ,  is  defined  by  the  expression: 

dW . 


T±  (1  +  n.¥7  (^)  Mi2) 


(A-2) 


where  n  is  1/2  for  laminar  flow  and  1/3  for  turbulent  flow.  Equating 
the  heat  flow  out  of  the  control  volume  for  a  steady  state  process  gives: 

(A-3) 

;v° 

or  i 

i  -i  ( A— 4 ) 

2 it  ri°  q„  +  2*  rj°  q  -  0 

where  jQ  is  either  0  or  1  for  two-dimensional  or  axi symmetric  flow, 
respectively. 
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Applying  Equation  A-l  to  the  values  of  the  individual  heat  fluxes  and 
solving  for  the  wall  temperature  gives: 

T  -  - 3 - - - -  . 

w  (1  +  Or./r^Nh./!!  )] 


For  flows  involving  a  constant  boundary  layer  edge  velocity  and  a 
constant  wall  temperature,  the  conductances,  hi ,  can  be  defined  by  the 
following  analytic  expressions: 


h±  -  0.332  Cp  p±  Re”J/2  Pr‘2/3 
for  laminar  flow,  and 

h±  -  0.0295  Cp  p±  u±  Pr_2/5 


(A-6) 


(A-  7 ) 


for  turbulent  flow  (Reference  40).  All  of  the  quantities  on  the  right- 
hand-side  of  Equation  A-5  are  then  known,  and  a  value  for  the  wall 
temperature  can  be  calculated  based  on  the  states  of  the  two  flowfields 
adjacent  to  the  nozzle  wall. 


AFWAL-TR-81 -31 61 


APPENDIX  B 

ADAPTIVE  FINITE  DIFFERENCE  MESH 

It  is  desirable  that  the  fine  mesh  region  of  the  computational  grid 
remain  in  the  areas  of  relatively  high  velocity  and  temperature  gradients 
as  the  solution  progresses  towards  convergence.  Hirt  (Reference  21) 
has  used  a  technique  in  the  solution  of  free  surface  flows  that  allows 
the  grid  to  adapt  as  the  solution  progresses.  The  following  kinematic 
equation  is  applied  in  the  region  where  the  nozzle  wake  and  shear  layer 
develop: 


f.cA(v-„f,  (B-U 

This  equation  ensures  the  condition  that  as  the  solution  converges, 
the  physical  slope  of  the  constant  n  finite  difference  cell  boundaries  is 
the  same  as  that  of  the  velocity  vectors  near  each  cell. 

Equation  B-l  can  be  converted  to  the  following  finite  difference 
form  for  application  to  a  computational  mesh: 


rM  ■  r".J  +  CA  “I’L 


n  n 

un  (AiJ  ~ 

i,j  n  n 

Xi,j  '  Xi-l,j 


)) 


(B-2) 


where 


.n 

Ui,j 


fi.j  *  1/4  (V 


n 

,  n 
Ui+l,j 

,  n  ,  n  . 

+  “i.j-i * 

(B-3) 

n 

.  n 

n  n 

i-l.J 

Vi+l,J 

vi,j-i +  ''i.j+P 

(B-4) 

and  where  is  a  constant  used  to  damp  the  grid  motion  with  respect  to 
time.  Spatial  averaging  of  the  velocity  components  is  applied  in  order 
to  reduce  the  effects  of  numerical  velocity  fluctuations  at  individual 
mesh  points  in  the  flow.  The  upwind  difference  form  of  the  cell  aspect 
ratio  term  (Ar/Ax)  is  also  utilized  to  achieve  better  stability  in  the 
finite  difference  equation. 
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Equation  B-l  is  applied  in  the  wake  region  of  the  flowfield  for  a 
line  of  constant  n  (j=constant) ,  where  the  specific  value  of  n  corresponds 
to  the  nozzle  inner  wall  for  overexpanded  flowfields,  and  to  the  nozzle 
outer  wall  for  underexpanded  flowfields.  Once  the  position  of  this  grid 
reference  line  is  established,  the  fine  mesh  region  corresponding  to  the 
wall  thickness  is  computed.  The  exponential  stretching  scheme  discussed 
in  Section  III  is  then  applied  for  each  value  of  £  in  the  regions  above 
and  below  the  fine  mesh  region  as  shown  in  Figure  41.  The  first  two  grid 
points  above  the  centerline  were  kept  fixed  at  constant  heights  for  all 
values  of  x.  This  prevented  large  numerical  errors  in  axisymmetric  cases 
involving  the  differencing  of  terms  containing  (1/r),  where  r  is  a  very 
small  number. 

The  constant  was  specified  in  the  range  of  0.3  -  0.6  in  order  to 
allow  the  grid  to  adapt  smoothly  as  the  solution  converged.  Larger  values 
of  caused  undesirable  oscillatory  motion  of  the  grid  reference  line 
with  respect  to  time. 

The  adaptive  grid  scheme  was  applied  once  during  every  iteration  of 
the  solution  algorithm.  The  number  of  points  allowed  to  "float"  on  the 
grid  reference  line  could  be  varied  during  the  course  of  the  solution. 

This  capability  was  utilized  primarily  during  the  start-up  portion  of  the 
numerical  solution,  where  only  a  limited  number  of  points  close  to  the 
nozzle  were  allowed  to  "float"  until  the  shear  layer  was  established. 

After  the  position  of  the  shear  layer  across  the  complete  computational 
domain  was  completely  established,  the  adaptive  grid  scheme  could  be 
turned  off  in  order  to  save  computer  time  during  the  remainder  of  the 
solution. 
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APPENDIX  C 

TWO-DIMENSIONAL  FLAT  PLATE  FAR  WAKE  SOLUTION 

In  order  to  test  the  validity  of  the  turbulence  model  in  the  far 

wake  region,  a  check  case  possessing  known  experimental  data  in  this 

region  was  computed.  The  experimental  data  for  this  case  is  that  of 

Toyoda  and  Hiriyama  (Reference  33),  which  involved  a  two-dimensional  thin 

flat  plate  at  a  Mach  number  of  1.60.  The  flat  plate  possessed  a  thickness 
-3  -4 

of  3.28x10  0  ft  (1  mm)  and  a  trailing  edge  thickness  of  3.28x10  ft 

(0.1  mm). 

The  following  flow  quantities  were  given  in  the  experimental  data: 

M  -  1.60 

QO 

P.  *  3.2  atm 

oo 

Re0  -  5300 

A  free  stream  stagnation  temperature  of  518. 7°R  was  assumed,  since  this 
value  was  not  given  experimentally.  The  wall  temperature  of  the  plate 
was  also  set  to  518. 7°R.  Using  the  previous  quantities,  the  boundary 
layer  momentum  thickness  at  the  trailing  edge  was  computed  as: 

0TE  -  2 . 95xlO~4  ft 

The  computational  solution  used  the  same  numerical  solver  (MacCormack's 
explicit  method)  as  that  of  the  nozzle  solutions.  The  eddy  viscosity 
models  utilized  were  identical  with  those  used  in  the  boundary  layer  and 
far  wake  regions  of  the  axisymmetric  nozzle.  The  boundary  conditions 
utilized  also  closely  resemble  those  of  the  nozzle.  The  inflow  conditions 
were  set  by  the  experimental  data  and  remained  fixed  for  the  duration  of 
the  solution.  The  upstream  u  component  of  velocity  near  the  trailing  edge 
of  the  plate  was  matched  to  that  of  the  experimental  data  as  shown  in 
Figure '42.  input  conditions  at  the  upstream  boundary  were  then  set  as: 

u(xlt  y,  t)  -  0.6295  uo<y/GTE)1/5  (C-l) 

v(x1,  y,  t)  -  0.0  (C-2) 


99 


AFWAL-TR-81-3161 


P(V  y’  tJ  “  P-/{1  +  (ar'>Mit1  -  Cu(xi,  y,  t)/uj2]} 
p(x  ,  y,  t)  =  p 

X  oo 


(C-3) 

(C-4) 


Both  the  upper  and  lower  freestream  boundaries  utilize  the  characteristic 
condition  applied  at  the  upper  boundary  of  the  nozzle  solutions.  The 
outflow  and  wall  boundaries  used  conditions  identical  to  those  utilized 
for  the  nozzle.  The  computational  solution  was  initially  started  by 
applying  the  upstream  profile  across  the  computational  domain: 


U(x,  y,  o)  =  U(x1§  y,  o)  (C-5) 

This  calculation  utilized  a  39x34  computational  mesh,  with 

exponential  mesh  stretching  employed  in  both  the  x  and  y  directions  as 

shown  in  Figure  43.  Minimum  grid  spacings  in  the  x  and  y  directions  were 
-4  -4 

4.00x10  ft  and  3.28x10  ft,  respectively.  The  physical  dimensions  of 
the  computational  flowfield  in  the  x  and  y  directions  were  0.125  ft  and 
0.070  ft,  respectively.  The  rate  of  data  processing  on  a  CDC  Cyber  175 
computer  was  0.0014  sec  per  grid  point  for  each  iterative  time  step. 

The  solution  was  computed  for  a  duration  of  four  characteristic  times 
(3200  iterations),  at  which  time  no  significant  change  was  detected  in 
the  dependent  variables.  The  result  was  then  taken  to  be  the  asymptotic 
solution. 

The  computational  solution  demonstrated  that  the  eddy  viscosity 
model  gives  good  results  in  the  far  wake  region.  The  maximum  velocity 
defect  generated  computationally  is  in  good  agreement  with  the  experi¬ 
mental  data  as  shown  in  Figure  21.  Figure  44  shows  that  the  velocity 
field  generated  computationally  evolves  from  the  boundary  layers  on 
the  plate  to  a  classic  wake  solution  very  rapidly  due  to  the  turbulent 
nature  of  the  flow. 
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e  43.  Computational  Mesh  Used  in  the  Flat  Plate  Solution 


Figure  44.  Computed  Flat  Plate  Velocity  Profiles 
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APPENDIX  D 

TWO-DIMENSIONAL  WEDGE-FLAT  PLATE  NEAR  WAKE  SOLUTION 

In  order  to  test  the  validity  of  the  turbulence  model  utilized  in 
the  near  wake  region  of  the  coflowing  nozzle  solutions,  a  check  case 
exhibiting  similar  physical  characteristics  and  possessing  known 
experimental  data  in  this  region  was  solved  numerically.  The  experimental 
case  selected  for  this  validation  study  is  one  obtained  by  Rom,  Seginer 
and  Kronzon  (Reference  35)  for  a  two-dimensional  wedge-flat  plate  in  a 
turbulent  supersonic  flowfield.  The  model  used  for  this  study  consisted 
of  a  sharp  15°  half  angle  wedge-flat  plate  with  a  base  height  of  10  mm 
and  a  chord  of  44  mm.  The  following  flow  conditions  were  given  in  the 
experimental  data: 

M  -  2.25 
00 

Re  -  1.5  x  106 
c 

PQ  -  40  psig 
00 

T  «  492°  R 
o 

oo 

A  computed  adiabatic  wall  temperature  of  467°R  based  on  the  flow  condition 
at  the  flat  plate  portion  of  the  model  was  used  in  the  computational 
solution. 

The  numerical  solution  used  the  same  computational  solver  and 
turbulence  model  as  that  of  the  coflowing  nozzle  solutions.  The  boundary 
and  initial  conditions  were  also  identical  to  those  used  in  the  coflowing 
nozzle  with  the  exception  of  the  jet  centerline  condition.  As  shown  in 
Figure  45,  this  condition  was  replaced  with  a  lower  freestream  boundary 
which  utilized  a  characteristic  scheme  similar  to  that  of  the  upper 
freestream  bounoary.  Since  a  value  for  the  boundary  layer  thickness  near 
the  trailing  edge  was  given  experimentally  ( 6 /h  =  0.15),  a  two-dimensional 
boundary  layer  code  was  used  to  generate  the  velocity  and  temperature 
profiles  on  the  upstream  boundary.  The  boundary  layer  starting  length 
was  adjusted  to  give  the  correct  boundary  layer  thickness.  One  additional 
condition  was  imposed  on  the  line  of  symmetry  for  the  wedge  (j  =  23). 


SIMPLE  WAVE  ASSUMPTION 

T  CONSTANT  ALONG  LINES  OF  CONSTANT  (0+u) 
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At  this  line  of  symmetry,  a  zero  v  velocity  component  was  enforced  to  help 
stabilize  the  wake  during  the  startup  of  the  numerical  solution  and  help 
accelerate  convergence  by  damping  any  numerical  shear  layer  oscillations 
in  the  wake. 

The  solution  was  calculated  using  a  45  x  45  point  computational  mesh 
with  exponential  stretching  employed  in  both  the  x  and  y  directions  as 
shown  in  Figure  46.  The  physical  dimensions  of  the  computational  flowfield 
in  the  x  and  y  directions  were  10cm  and  7cm,  respectively.  Minimum  grid 
spacing  in  the  x  and  y  directions  was  set  at  0.5  mm.  This  gave  a  value  of 
u/ug  =  0.83  for  the  first  point  in  the  boundary  layer  above  the  nozzle 
wall,  which  corresponds  to  an  identical  value  in  the  jet  boundary  layer 
of  the  axisymmetric  nozzle  solutions.  Thus,  truncation  error  should  be 
similar  in  this  region  for  both  the  wedge-flat  plate  and  the  nozzle.  The 
rate  of  data  processing  on  a  CDC  Cyber  175  computer  was  0.0014  sec  per 
grid  point  for  each  iterative  time  step.  The  solution  w^s  allowed  to 
progress  for  approximately  four  characteristic  timesteps  at  which  time 
the  change  in  the  dependent  variables  was  less  than  0.5 %  per  characteristic 
time  period.  This  condition  was  then  considered  to  be  the  converged 
asymptotic  solution. 

Specific  features  of  the  experimental  flowfield  in  the  near  wake  were 
reproduced  in  the  computational  solution  and  can  be  distinguished  in  the 
plots  of  Mach  number  contours  and  velocity  profiles  shown  in  Figures  47 
and  48.  Several  flow  features  which  were  numerically  observed  include  the 
existence  of  the  boundary  layers  along  the  horizontal  walls  of  the  body, 
the  turning  of  the  flow  through  the  corner  expansion  fans,  the  existence 
of  the  subsonic  recirculating  "deadwater"  region  adjacent  to  the  base  of 
the  body,  flow  recompression  through  the  trailing  shocks,  and  the  evolution 
of  the  wake  to  a  classic  far  wake  flow.  The  weak  lip  shock  evident  in 
the  experimental  data  was  not  readily  evident  in  the  numerical  solution. 
This  may  be  attributed  to  the  fact  that  the  numerical  method  tends  to 
smear  shocks,  and  thus  has  difficulty  locating  shocks  which  are  very  weak. 

Quantitative  accuracy  of  the  numerical  solution  is  identified 
through  the  use  of  the  given  experimental  static  pressure  and  pitot 
pressure  data.  A  comparison  of  the  axial  static  pressure  distribution 
along  the  line  of  symmetry  is  shown  in  Figure  49.  The  computed  static 
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Figure  46.  Computational  Mesh  Used  for 
Wedge-Flat  Plate  Solution 
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pressure  distribution  is  within  5%  of  the  measured  values  except  in  the 
region  of  reci rculating  flow  (x  £  1.0cm),  where  there  is  up  to  a  7% 
discrepancy  between  the  experimental  and  computed  values.  However, 
static  pressure  probes  like  the  one  used  to  obtain  the  data  are  very 
sensitive  to  flow  angularity,  and  thus  less  reliable  in  areas  of 
recirculating  flow.  A  more  reliable  comparison  is  that  of  the  pitot- 
pressure  surveys  shown  in  Figure  50.  These  measurements  are  less 
susceptible  to  errors  resulting  from  flow  angularity.  This  figure  shows 
excellent  agreement  in  the  "deadwater"  region  of  recirculating  flow. 

The  numerical  solution  smears  the  weak  beginning  of  the  recompression 
trailing  shocks  at  0.5cm,  but  correctly  simulates  them  at  the  correct 
values  further  downstream  (x  >_  1.0).  This  figure  particularly  demonstrates 
that  the  phenomena  present  in  the  near  wake  are  accurately  simulated  by  the 
present  computational  method  and  turbulence  modeling. 


Figure  50.  Pitot  Pressure  Profiles  in  the  Near  Wake  of  the  Two- 

Dimensional  Wedge-Flat  Plate  (Symbol s-Experimental  Data 
(Reference  34),  Solid  Lines  -  Computational  Solution) 
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APPENDIX  E 

AXISYMMETRIC  NOZZLE  SOLUTION  SIMULATING  INTERNAL 
BOUNDARY  LAYER  SEPARATION 

The  peculiar  reversal  in  the  nozzle  base  pressure  coefficient  as 
the  pressure  ratio  was  lowered  below  P./P^  =  0.18  was  believed  to  be 
caused  by  boundary  layer  separation  in  the  divergent  portion  of  the  nozzle 
(Figure  38).  This  was  further  substantiated  through  an  examination  of  the 
schlieren  photographs  (Reference  16)  which  show  a  definite  change  in  flow 
character  near  this  value  (Figure  39).  To  test  this  hypothesis  a  numerical 
solution  was  computed  for  a  case  where  the  jet  total  pressure  corresponded 
to  that  of  the  attached  strong  shock  solution  (P^/P*,  =  0.15),  but  where 
the  jet  Mach  number  equals  2.60  based  upon  an  isentropic  expansion  in  the 
nozzle.  For  this  jet  Mach  number  a  reduction  in  nozzle  area  ratio  was 
assumed  from  A/A*  =  4.2  at  Mach  3  to  A/A*  =  2.9  at  Mach  2.60.  This 
reduced  expansion  rate  was  intended  to  roughly  simulate  the  effects  of 
boundary  layer  separation  in  the  nozzle,  while  retaining  grid  geometry 
and  fineness  identical  to  the  previously  calculated  attached  jet  flow 
case. 

The  following  jet  flow  parameters  were  used  in  this  solution: 

M,  *=  2.60  T  =  551 .0°K 

J  w 

T  .  -  580. 5°R 

P  .  =  1636  psf 

oj  H 

These  jet  conditions  produce  an  actual  nozzle  pressure  ratio  of 

P./P  =  0.276  versus  the  calculated  value  of  P./P  =  0.150  which  was 
J  00  j  °° 

assumed  in  Reference  16. 

This  solution  was  initialized  using  the  solution  for  the  attached 
strong  shock  case  with  a  calculated  nozzle  pressure  ratio  equal  to  0.150. 
Except  for  the  upstream  jet  boundary,  boundary  conditions  and  mesh 
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configuration  identical  to  those  of  the  attached  case  were  utilized. 

At  the  upstream  jet  boundary  the  following  velocity  profile  was  fixed: 

u  -  2002  ft/sec  0  <  r/r,  _  <  0.835 

u  -  0  0.835  <  r/rjet  <  1 

v  -  0  0  <  r/r.  <  1 

—  jet  — 

As  shown  in  Figure  51,  this  velocity  input  profile  increased  the 
displacement  thickness  of  the  jet  boundary  layer  in  rough  approximation 
to  that  caused  by  separation  in  the  nozzle.  In  this  case  the  jet  flow 
is  assumed  to  be  turbulent  in  nature. 

As  shown  in  Figure  52,  a  strong  shock  solution  was  obtained,  with 
a  reflection  length  of  *s/rjet  =  1-29,  about  2%  greater  than  the  experi¬ 
mental  value  at  this  jet  total  pressure  case.  The  base  pressure 
coefficient  obtained  was  equal  to  -0.264,  about  8 %  less  than  the  experi¬ 
mental  value  and  in  much  better  agreement  than  that  of  the  attached  jet 
boundary  layer  case. 

Although  discrepancies  in  this  solution  such  as  the  smaller  diameter 
of  the  computational  Mach  disc  and  the  appearance  of  a  physically  non- 
existant  secondary  Mach  disc  further  downstream  are  apparent,  this  case 
does  demonstrate  that  separation  and  its  effect  on  the  flow  expansion 
in  the  nozzle  can  significantly  impact  the  resultant  base  pressure 
coefficient  of  the  nozzle.  Therefore,  future  computations  for  low  pres¬ 
sure  ratios  should  commence  at  the  throat  section  to  insure  that  nozzle 
separation  is  properly  considered. 


Ill 


r  Contours  for  the  Separated  FI 
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APPENDIX  F 

INVESTIGATION  OF  NUMERICAL  ERROR 

Two  error  analyses  were  performed  on  the  numerical  technique  used 
to  obtain  solutions  for  the  nozzle  flowfields.  The  first  consisted  of 
examining  the  error  generated  by  the  numerical  algorithm  used  to  solve  the 
Navier-Stokes  equations.  The  second  was  a  study  of  the  effect  of  re¬ 
positioning  first  the  downstream  boundary,  and  subsequently  the  upper 
boundary  to  regions  containing  only  minor  flow  gradients  normal  to  each 
boundary.  These  analyses  are  discussed  in  detail  in  the  following  sections. 

1.  TRUNCATION  ERROR  ANALYSIS 

MacCormack's  explicit  finite  difference  algorithm  is  an  equivalent 
second  order  accurate  numerical  solver.  The  final  converged  solution  for 
any  case  computed  by  this  algorithm  should  then  satisfy  the  Navier-Stokes 
equations  at  all  interior  node  points  with  second  order  accuracy.  A 
numerical  check  on  this  accuracy  was  conducted  using  a  typical  converged 
solution  and  the  following  procedure. 

As  shown  in  Section  II,  the  axisymmetric  Navier-Stokes  equations 
can  be  written  as: 


j>U  .  3F  1  3(rG)  H 
3t  3x  r  3r  r 


(F-l) 


A  nondimensional ized  finite  difference  formulation  of  these  equations 
can  then  be  written  as: 
t  ...  t 


(_£)  iLL  +  r^F  1  MrG)  H.  _ 

Si  }  At  +  U  [Ax  +  7  ~JT  ~  P  "  E  “  Error 


(F-2) 


where  U  is  defined  as: 

oo 
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and  t  =  5  x  10"^  sec  is  set  to  insure  that  the  lead  terms  on  the  left 

c 

hand  side  of  Equation  F-2  are  of  order  one. 

The  particular  check  case  selected  was  that  containing  the  strong 
Mach  disc  shock  structure  (P./P^  =  0.150),  since  this  case  contains  both 
substantial  regions  of  subsonic  and  supersonic  flow,  and  oblique  as  well 
as  normal  shock  waves.  The  MacCormack  solution  for  this  case  was  used  as 
input  for  Equation  F-2  where  the  left  hand  side  of  the  equation  was 
computed  at  all  interior  grid  points  using  a  standard  two-dimensional 
second  order  central  differencing  scheme  applied  on  the  transformed 
computational  plane.  The  magnitude  of  the  Error  vector  (E)  is  then  an 
indication  of  how  close  MacCormack' s  method  is  to  an  alternate  second 
order  accurate  solver.  The  following  root  mean  square  (RMS)  values  of  E 
were  obtained  over  the  interior  of  the  computational  domain: 


’  1 

E 

p 

0.026' 

E 

K 

0.020 

pu 

E 

0.011 

pv 

E 

0.026 

pe 

This  result  indicates  that  over  the  domain,  MacCormack's  algorithm  and 
the  two-dimensional  central  difference  scheme  are  equivalent  to  within 
three  percent. 

Finally,  the  Error  vector  E  was  examined  over  the  computational  domain 
to  determine  which  regions  generated  the  highest  magnitude  of  error. 

Since  the  RMS  error  value  of  the  continuity  equation  was  one  of  the 
largest,  and  the  error  distribution  was  representative  of  that  in  the 
other  equations,  it  is  shown  in  Figure  53.  In  this  figure,  only  error 
values  greater  than  the  RMS  value  are  shown  as  contours.  Regions  con¬ 
taining  the  largest  error  consist  of  those  containing  shock  waves  and 
that  containing  the  expansion  fan  near  the  sharp  corner  of  the  nozzle. 

In  these  shock  regions,  strong  flow  gradients  exist  over  areas  with 
fairly  coarse  finite  difference  mesh  spacing.  Although  the  wake  region 
also  contains  strong  gradients  within  the  mixing  layer,  it  lies  within 
a  fine  mesh  region  of  the  grid  that  produces  much  less  numerical  error. 
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RMS  VALUE  (0.03  MIN,  0.12  MAX,  0.03 
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This  analysis  further  demonstrates  the  desirability  of  utilizing 
adaptive  mesh  schemes  that  can  align  the  grid  with  flow  gradients  as  the 
solution  progresses  to  convergence. 

2.  BOUNDARY  POSITION  ANALYSIS 

Although  the  upstream  boundary,  the  centerline  boundary,  and  the 
position  of  the  nozzle  walls  were  fixed  by  the  definition  of  the  problem 
to  be  solved,  the  placement  of  both  the  downstream  boundary  and  the  upper 
boundary  was  left  to  the  discretion  of  the  computational  investigator. 

It  was  desirable  to  place  these  boundaries  as  close  to  the  nozzle  as 
possible  in  order  to  achieve  computational  efficiency  in  a  compact 
domain.  In  the  axisymmetric  nozzle  solutions,  both  of  these  boundaries 
were  located  in  regions  in  which  flow  gradients  existed  due  to  the 
presence  of  shock  and  expansion  waves  and  viscous  phenomena  such  as 
shear  layers  and  wakes  in  the  flowfield.  The  assumption  is  made  that 
positioning  these  two  boundaries  in  flow  gradient  regions  does  not  affect 
the  computational  solutions  obtained.  To  validate  this  assumption,  each  of 
these  boundaries  was  repositioned  a  greater  distance  from  the  nozzle  to 
regions  containing  only  minor  flow  gradients  normal  to  each  boundary.  The 
resulting  effect  on  the  shock  wave  structure  as  well  as  on  the  nozzle 
base  pressure  coefficient  in  the  numerical  solution  was  then  observed. 


The  nozzle  case  at  which  P./P  =  0.251  was  examined  for  this 

J  00 

particular  study.  The  downstream  boundary  contains  primarily  supersonic 
outflow  with  an  embedded  wake  region  of  subsonic  outflow.  The  downstream 
boundary  was  extended  from  its  original  position  at  x/rjet  =  3.0  to  a 
new  value  of  x/rje<-  =  6.0.  This  stretching  was  achieved  by  the  addition 
of  twelve  grid  points  axially  to  the  original  mesh.  The  outflow  at  this 
new  position  was  totally  supersonic  in  nature  with  only  minor  gradients 
in  existence  normal  to  the  boundary.  As  shown  in  Figure  54,  no  changes 
were  evident  in  the  shock  structure  contained  in  the  original  domain. 

The  computational  base  pressure  coefficient  remained  unchanged  at  a 
value  of  Pg  =  -0.299. 
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The  effect  of  repositioning  the  upper  boundary  to  a  radial  distance 
at  which  flow  gradients  are  not  present  was  then  examined.  The  upper 
boundary  of  the  axially  stretched  case  was  extended  from  its  original 


value  of  r/r,ej.  =  3.0  to  a  new  value  of  r/rje^.  =  6-0  as  shown  in 
Figure  55.  In  this  case  eight  grid  points  were  added  radially  to  the 


axially  stretched  mesh.  Again,  no  changes  were  evident  in  the  numerical 


shock  structure,  and  the  computational  value  of  the  nozzle  base  pressure 


coefficient  remained  at  Pg  =  -0.299. 


Since  repositioning  these  boundaries  had  essentially  no  effect  on 
the  computational  solution  that  was  examined,  the  application  of  these 
boundary  conditions  in  the  original  regions  containing  mixed  supersonic- 
subsonic  flow  on  the  outflow  boundary  and  substantial  flow  gradients 
on  both  bound  ries  was  valid. 
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